A  new  certification  framework  for  the  port  reduced  static  condensation 

reduced  basis  element  method 


Kathrin  Smetana 

Department  of  Mechanical  Engineering,  Massachusetts  Institute  of  Technology,  77  Massachusetts  Avenue,  Cambridge, 

M A- 02 139,  USA 


Abstract 

In  this  paper  we  introduce  a  new  certification  framework  for  the  port-reduced  static  condensation  reduced 
basis  element  (PR-SCRBE)  method,  which  has  been  developed  for  the  simulation  of  large  component  based 
applications  such  as  bridges  or  acoustic  waveguides.  In  an  offline  computational  stage  we  construct  a  library  of 
interoperable  parametrized  reference  components;  in  the  subsequent  online  stage  we  instantiate  and  connect 
the  components  at  the  interfaces/ports  to  form  a  system  of  components.  To  compute  a  “truth”  finite  element 
approximation  of  the  (say)  coercive  elliptic  partial  differential  equation  on  the  component  based  system  we 
use  a  domain  decomposition  approach.  For  an  efficient  simulation  we  employ  two  different  types  of  model 
reduction  —  a  reduced  basis  (RB)  approximation  within  the  interior  of  the  component  (Huynh  et  al.  2013) 
and  empirical  port  reduction  (Eftang  and  Patera  2013)  on  the  ports  where  the  components  connect. 

We  demonstrate  the  well-posedness  of  the  PR-SCRBE  approximation  and  introduce  a  new  certification 
framework.  To  assess  the  quality  of  the  port  reduction  we  use  conservative  fluxes.  We  adapt  the  standard 
estimators  from  RB  methods  to  the  SCRBE  setting  to  derive  an  a  posteriori  error  estimator  for  the  RB- 
error  contribution.  In  order  to  combine  the  a  posteriori  estimators  for  both  error  contributions  and  derive 
a  rigorous  a  posteriori  error  estimator  for  PR-SCRBE  we  adapt  techniques  from  multi-scale  methods  and 
component  mode  synthesis.  Finally,  we  prove  that  the  effectivity  of  the  derived  estimator  can  be  bounded. 
We  provide  numerical  experiments  for  heat  conduction  and  linear  elasticity  to  show  that  the  derived  a 
posteriori  error  estimator  provides  an  effective  estimator.  Moreover  we  demonstrate  the  applicability  of  the 
introduced  certification  framework  by  analyzing  the  computational  (online)  costs. 
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1.  Introduction 

Within  many  engineering  applications  the  considered  structure  allows  for  a  natural  decomposition  in  com¬ 
ponents.  Examples  are  bridges,  buildings,  aircrafts,  oil  and  gas  platforms,  musical  instruments  or  mufflers. 
One  popular  method  for  the  simulation  and  analysis  of  such  large  engineering  systems  is  component  mode 
synthesis  (CMS).  The  CMS  approach  introduced  in  [1,  2]  uses  the  eigenmodes  of  local  constrained  eigenvalue 
problems  for  the  approximation  within  the  interior  of  the  component  and  static  condensation  to  arrive  at  a 
(Schur  complement)  system  associated  with  the  coupling  modes  on  the  interfaces  or  ports.  In  more  recent 
works  also  the  static  or  coupling  modes  are  chosen  as  eigenmodes  [3,  4,  5]  and  used  within  an  adaptive  scheme 
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based  on  a  posteriori  error  estimators  [5] .  One  drawback  of  the  CMS  approach  is  the  rather  slow  convergence 
of  eigenmodal  expansions.  In  contrast  the  reduced  basis  element  (RBE)  method  [6]  employs  a  reduced  basis 
expansion  [7]  within  each  component  or  subdomain  and  Lagrange  multipliers  to  couple  the  local  bases  and 
hence  compute  a  global  solution  of  the  considered  parameter  dependent  partial  differential  equation  for  each 
admissible  parameter.  The  RBE  method  thus  profits  from  the  fact  that  RB  approximations  yield  a  rapid 
and  in  many  cases  exponential  convergence  [8] . 

A  combination  of  RB  methods  and  domain  decomposition  approaches  has  for  instance  also  been  considered 
in  [9,  10].  Similarly  RB  methods  have  been  employed  in  the  framework  of  a  multi-scale  finite  element  method 
to  construct  local  reduced  spaces  for  the  approximation  of  fine-scale  features  on  the  coarse  grid  elements  in 
[11,  12,  13],  where  the  latter  corresponds  to  the  "components"  in  the  RBE  method. 

In  this  paper  we  consider  the  port  reduced  static  condensation  reduced  basis  element  (PR-SCRBE)  method 
[14,  15]  which  combines  ideas  of  the  CMS  and  the  RBE  method.  SCRBE  as  introduced  in  [14]  has  been 
successfully  applied  amongst  others  to  heat  exchanger  systems  [16]  and  musical  instruments  and  mufflers 
[17]  and  PR-SCRBE  introduced  in  [15]  to  heat  transfer  [15]  and  large  structures  in  linear  elasticity  [18]. 
A  key  advantage  of  the  PR-SCRBE  method  is  its  "bottom-up"  approach  from  an  offline-prepared  library 
of  completely  interchangeable  and  interoperable  parametrized  components  to  many  different  online-formed 
global  systems  of  components.  To  construct  rapid  convergent  reduced  bases  we  use  an  RB  approximation  in 
the  component  interior  [14]  and  empirical  port  reduction  based  on  RB  techniques  [15]  on  the  ports,  while  the 
coupling  is  realized  by  applying  static  condensation.  More  precisely,  we  first  compute  coupling  modes  as  the 
harmonic  liftings  of  the  empirical  port  modes.  Subsequently,  we  compute  RB  approximations  of  the  compo¬ 
nent  interior  "bubble"  functions,  each  associated  with  a  different  coupling  mode,  which  solve  the  parameter 
dependent  PDE  within  each  component  and  hence  account  for  material  or  geometric  parameters.  Employing 
a  different  RB  approximation  for  each  bubble  function  yields  a  rapid  convergent  and  computationally  efficient 
approximation . 

The  main  contribution  of  this  paper  is  the  introduction  of  a  new  certification  framework  for  the  PR-SCRBE 
method  which  provides  a  rigorous  and  effective  a  posteriori  error  estimator  of  the  error  between  the  PR- 
SCRBE  approximation  and  the  "truth"  finite  element  approximation  on  the  global  system  in  an  energy 
norm.  We  adapt  techniques  from  RB  methods  [7]  to  derive  an  a  posteriori  error  estimator  for  the  RB  error  of 
the  bubble  approximation  in  the  component  interior  by  considering  weighted  Riesz  representations.  To  derive 
an  a  posteriori  error  estimator  for  port  reduction  we  adapt  the  concept  of  conservative  fluxes  introduced  in 
Hughes  et  al.  [19]  to  the  PR-SCRBE  setting.  We  compute  weak  fluxes  of  the  PR-SCRBE  approximation  on 
each  port  with  respect  to  the  full  port  space  and  employ  the  jumps  to  define  an  a  posteriori  error  estimator 
for  port  reduction.  Thanks  to  the  weak  flux  continuity  of  the  "truth"  solution  the  proposed  estimator  thus 
measures  how  well  the  PR-SCRBE  approximation  behaves  like  the  "truth"  solution  on  the  ports.  We  remark 
that  the  concept  of  conservative  fluxes  has  been  used  as  an  error  indicator  in  the  formulation  of  an  adaptive 
variational  multi-scale  method  in  [20]  and  an  adaptive  multi-scale  finite  element  method  in  [21]  to  determine 
the  computational  domains  for  the  local  fine  scale  problems.  To  combine  both  error  estimators  we  exploit 
that  due  to  the  fact  that  the  coupling  modes  are  chosen  as  harmonic  extensions  the  span  of  the  coupling 
modes  and  the  span  of  the  bubble  functions  for  each  component  are  orthogonal.  Orthogonality  with  respect 
to  the  bilinear  form  has  been  exploited  for  the  derivation  of  an  a  posteriori  error  estimator  for  the  CMS 
approach  in  [5].  In  the  recently  introduced  local  orthogonal  decomposition  method  [22]  the  coarse  and  fine 
scale  spaces  are  constructed  in  such  a  way  that  orthogonality  with  respect  to  the  bilinear  form  is  obtained 
and  then  exploited  for  an  error  analysis. 

Our  a  posteriori  error  estimator  features  several  important  innovations.  First,  in  contrast  to  earlier  con¬ 
tributions  [14,  15,  18],  it  provides  an  estimator  for  the  energy  norm  of  the  error.  Second,  its  validity  does 
not  rely  on  the  quality  of  the  RB  approximation,  whereas  the  estimator  proposed  in  [14]  provides  only  a 
rigorous  bound  if  the  sum  over  the  RB-error  in  the  entries  of  the  Schur  complement  matrix  is  smaller  than 
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the  smallest  eigenvalue  of  the  “truth”  Schur  complement  matrix.  Moreover,  the  error  estimator  is  based 
on  local  component-wise  element  indicators  and  may  hence  be  employed  within  an  adaptive  PR-SCRBE 
scheme.  Note  that  in  contrast  the  a  posteriori  error  estimator  proposed  for  the  CMS  approach  in  [5]  relies 
on  a  weak  residual  with  respect  to  the  global  port  space.  For  the  localized  reduced  basis  multi-scale  method 
an  a  posteriori  bound  using  local  error  indicators  based  on  conforming  flux  reconstruction  is  stated  in  [23] . 
We  also  derive  upper  bounds  for  the  effectivities  both  of  the  error  estimator  and  the  local  indicators.  Finally, 
we  demonstrate  in  numerical  experiments  that  the  derived  estimator  is  also  computationally  efficient  as  its 
computation  requires  only  relatively  small  additional  online  costs  in  comparison  to  the  computation  of  the 
PR-SCRBE  approximation. 

The  remaining  part  of  this  paper  is  organized  as  follows.  In  Section  2  we  describe  the  parametrized  com¬ 
ponent  based  static  condensation  framework  as  introduced  in  [14].  In  the  following  Section  3  we  present 
the  PR-SCRBE  method  [14,  15]  and  prove  well-posedness  of  the  PR-SCRBE  approximation.  The  main  new 
contribution  of  this  paper  is  developed  in  Section  4  where  we  derive  the  a  posteriori  error  estimator  and  prove 
its  effectivity.  Subsequently  we  discuss  the  computational  costs  both  of  the  PR-SCRBE  approximation  and 
the  introduced  error  estimator  in  Section  5.  Finally,  we  present  numerical  experiments  for  heat  conduction 
and  linear  elasticity  in  Section  6  to  validate  the  theoretical  results  and  draw  some  conclusions  in  Section  7. 

2.  Component  based  static  condensation 

The  key  concepts  of  the  PR-SCRBE  approximation  are  a  library  of  parametrized  archetype  components 
associated  with  reference  domains  and  a  parametrized  PDE  (§2.1.1)  and  a  system  of  instantiated  components, 
mapped  to  their  physical  coordinates  and  connected  at  ports  (§2.1.2).  Each  archetype  component  represents 
a  specific  geometric  form,  for  instance  a  beam  or  a  fin,  and  may  feature  various  ports  of  different  types.  To 
simplify  the  notation  we  will  however  consider  in  this  paper  only  one  archetype  component  and  one  port 
type. 

To  compute  a  numerical  approximation  of  the  parametrized  PDE  for  the  global  system  we  employ  a  (multi- 
domain)  domain  decomposition  approach  and  eliminate  degrees  of  freedom  in  the  component  interior  by 
static  condensation  (§2.2). 

2.1.  Component  library 

2.1.1.  Archetype  (reference)  component 

The  considered  archetype  component  is  associated  with  a  bounded  reference  domain  C  d  =  1,2,3 
with  Lipschitz  boundary  Oil.  It  is  equipped  with  P7  local  ports  jj  C  912,  j  =  1, ...,  P7  at  which  instantiations 
of  this  archetype  components  can  be  connected  in  the  online  stage  as  depicted  in  Fig.  1.  We  assume  that 
local  ports  are  mutually  separated  by  manifolds  E  C  dfl  with  positive  Hausdorff  measure.  For  the  treatment 
of  systems  of  components  with  intersecting  ports  we  refer  to  [4] . 

Furthermore,  we  associate  with  the  archetype  component  a  parametrized  PDE  modeling  the  physical  phenom¬ 
ena  of  the  target  application.  To  formulate  the  corresponding  variational  problem  we  introduce  an  archetype 
continuous  bilinear  form  a(-,-;/t)  :  X  x  X  — >  K.  and  a  continuous  linear  functional  /(•;/),)  :  X  — >  R,  where 
(HQ(fl))dr  C  X  C  (P1(12))dT'.  Here,  /t  £  V  C  Rz  denotes  a  parameter  vector  which  describes  for  instance  a 
material  parameter  as  Young’s  modulus  or  modifies  the  geometry  of  the  component  say  via  a  dilation  and 
T>  is  the  corresponding  2-dimensional  parameter  space.  Moreover,  dr  equals  one  for  scalar  valued  and  d  for 
vector  valued  parametrized  PDEs  and  Dirichlet  boundary  conditions  may  be  only  prescribed  on  ports.  We 
assume  an  affine  parameter  dependence  of  a(-,  •;  (i)  and  /(•;  fi),  i.e. 

Qa  Qf 

£(•>•;/})  =  J2^q(h)aq{-,-)  and  /(•;  p)  =  ^  (1) 

g=l  9=1 
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(a)  archetype  (reference)  component 


(b)  system  of  two  rotated  and  deformed  components 


Figure  1:  A  library  with  one  archetype  component  (a):  The  archetype  reference  domain  is  and  the  archetype  component  has 
P'y  =  3  local  ports  7^,  j  =  1, 3.  A  system  of  I  =  2  components  1  and  2  which  are  instantiations  of  the  archetype  (b):  The 
component  domains  are  £1  =  7i(f2)  and  £2  2  =  T2  (£2)  and  the  components  are  connected  at  the  global  port  T3.  Thus  we  have 
for  instance  tt\  =  {(1,1)},  7T2  =  {(l,2)},7r3  =  {(1,  3),  (2,  3)}  and  conversely  IIi(l)  =  1,  Hi (2)  =  2  and  IIi(3)  =  112(3)  =  3.  Both 
mappings  7i,  i  =  1,2  are  compositions  of  a  rotation  7}ro*  and  a  dilation  7^e^,  i  =  1,2  indicated  by  the  shading  in  gray.  Note 
that  thanks  to  the  application  of  'J™av  to  the  dependent  variables  we  obtain  a  field  in  the  physical  system  coordinates  which 
has  the  same  spatial  orientation  as  in  the  coordinate  system  of  the  archetype  reference  component. 


where  a9(-,  •)  :  X  x  A'  — >■  R,  q  =  1, Qa  and  fq(-)  :  X  — R,  q  =  1, Q?  are  parameter  independent  bilinear 
and  linear  forms,  and  0“  and  9 f  are  parameter  dependent  functions.  For  general  parameter  dependent 
forms  an  (empirical)  interpolation  [24]  may  be  invoked  to  realize  an  (approximate)  affine  decomposition. 
Furthermore,  we  define  an  inner  product  and  the  corresponding,  induced  norm  as 

(■>  Ox  ;  X  x  X  — ►  R  and  ||  •  ||^-  :=  \J (•,  Ox- 

We  require  ||  •  ||^  to  be  a  semi-norm  on  (H1(£l))dr  and  a  full  norm  on  the  space  {  v  £  (H1(Q))dr  :  v\^.  = 
0  for  at  least  one  j  S  {1,  ...,P7}}.  For  notational  convenience  we  also  introduce  a  (geometric)  parameter- 
dependent  inner  product  on  the  archetype  component  as  (■,•;  (i)  :  X  x  X  — >  R  and  the  induced  norm 

IMIx'/i  :=  ((^i  A)  ^or  ^  G  X.  Finally,  we  introduce  "truth"  archetype  finite  element  (FE)  spaces 

X-^r  C  X  of  dimension  J\f.  For  each  local  port  within  an  archetype  component  we  introduce  the  trace  space 
A  j  =  {v  £  H ■  v  =  v\xj.,  v  £  A'}  and  its  discrete  counterpart  —  the  discrete  port  space 

Afp  :=  XM\*tj  =span{xyi,...,xiiAA-p}>  i  =  l,...,P7,  (2) 

of  dimension  Af-p  with  basis  {Xj,k}'kZi- 
2.1.2.  System 

In  the  online  stage  we  form  the  target  structure,  say  a  bridge  or  an  electronic  module,  and  hence  a  global 
system  of  I  components  mapped  to  physical  (system)  coordinates.  Each  component  i,  i  =  1, ...,/  is  instan¬ 
tiated  from  one  archetype  component  as  assumed  in  the  beginning  of  this  section.  We  also  introduce  the 
parameter  vector  fit  £  T>i  C  V.  Note  that  T>i  differs  between  distinct  components  in  the  global  system  if  we 
prescribe  for  instance  a  traction  on  one  end  of  a  beam  structure  as  in  Subsection  6.2.  We  may  then  define 
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an  affine  parameter  dependent  geometric  transformation  7 $  :  from  (archetype)  reference  to  physical 

(system)  coordinates,  where  we  denote  with  f 2,  =  7 )(fi)  the  instantiated  and  transformed  component  domain. 
Accordingly  7 =  Tii'jj),  j  =  1,  ...,P7  are  the  transformed  ports.  We  assume  that  7)  is  a  composition  of  a 
deformation  7*^  and  a  rotation  7)rot,  i.e.  7)  =  T{otT^e^  with  {x)  =  Gi.x  + 1»,  x  £  fl.  Here,  Gj  G  Rdxd 
is  the  Jacobian  of  Tfe^  and  t,:  G  is  a  translation  vector.  Hence,  7)de^  may  describe  a  composition  of  (say) 
a  dilation  and  a  translation,  but  rotations  are  not  permitted.  In  this  paper  we  assume  that  flj  and  hence  71 
do  not  degenerate.  Moreover,  we  suppose  that  7)  when  applied  to  a  port  is  a  pure  rigid  body  transformation, 
i.e.  a  composition  of  a  rotation  and  translation,  which  is  for  instance  fulfilled  for  the  mappings  in  Fig.  1. 
Details  on  the  difficulties  that  arise  when  this  assumption  is  dropped  and  possible  solutions  are  provided 
below. 

Next  we  define  the  global  (system)  domain  fi  as  Cl  =  U,(=1  Q,:  and  recall  that  the  instantiated  and  trans¬ 
formed  components  are  connected  at  their  ports.  More  precisely  we  say  that  two  components  i  and  i'  are 
connected  at  a  global  port  Tp  if  n  IV  =  Tp.  We  register  those  connections  in  a  global-to-local  index  set 
%  =  {(*,  j),  which  contains  also  the  indices  of  the  respective  local  ports  7 ij  and  7 pji  at  which  the 

components  i  and  i'  connect.  Note  that  a  global  port  may  also  lie  on  dil  and  is  then  only  associated  with 
one  local  port  7 of  component  i:  hence  we  set  irp  =  {(*,_))}  in  this  case.  Additionally  we  introduce  for 
each  instantiated  component  i,  i  =  1, ...,/,  a  local-to-global  port  map  n^,  which  is  defined  as  n, (j)  =  p  if 
(i.  j)  G  7TP,  j  =  1  ,  ...,P7  and  p  =  1,...,P(^  as  illustrated  in  Fig.  1.  Here,  Pq  is  the  number  of  global  ports 
in  the  system.  Finally,  we  denote  with  Pr  <  Pq  the  number  of  global  ports  on  which  we  do  not  impose 
Dirichlet  conditions. 

To  define  a  "truth"  FE  solution  for  the  global  system  we  first  introduce  for  each  component  i,  i  =  1, ...,/  a 
mapped  FE  space 

Xj^  :=  span{7”map(u  o  7)-1),  v  £  XM},  (3) 

where  7)mop  =  T[ot  for  vector- valued  functions  and  T^nap  equals  the  identity  map  for  scalar- valued  functions. 
We  emphasize  that  we  should  instead  use  the  notation  X1^  (/q),  but  to  improve  readability  we  omit  here  and 
henceforth  the  argument  p  for  spaces  that  depend  (merely)  on  the  geometrical  parameters.  Next,  we  introduce 
for  i  =  1, ...,  I  and  j  =  1, ...,  P7  a  mapped  discrete  port  space 

AujT  ;=span  {Xi,j,k,  k  =  l,...,Afp}  with  Xi,j,k  ■=  0  T~l).  (4) 

To  obtain  continuity  of  the  global  solution  it  is  essential  to  require  compatible  port  spaces  in  the  sense  that 

Xi,j,k  =  Xi',j',k  for  7 rp  =  {(i,j),  («',/)}  (5) 

for  all  global  ports  rp  within  the  global  system  at  which  two  components  connect.  One  key  point  to  satisfy 
(5)  is  to  define  the  (archetype)  port  spaces  KJ^'P  according  to  the  orientation  of  the  respective  (archetype) 
port  7 j  within  the  component.  Moreover,  for  vector-valued  functions  the  application  of  7”rot  to  the  dependent 
variables  (see  (3)  and  (4))  ensures  the  right  orientation  of  the  mapped  port  space  in  the  global  system. 
We  refer  to  [18]  for  a  detailed  discussion  of  the  port  compatibility  requirement  (5)  in  the  case  of  vector-valued 
functions. 

We  may  now  introduce  a  global  solution  space  (Tf^fl))^  C  X  C  ( H1(Q))dr  and  a  global  FE  space  X1^  c  X 
as  X^  =  ® (=1  Xf .  We  only  consider  homogeneous  Dirichlet  boundary  conditions  as  non-homogeneous 
Dirichlet  boundary  conditions  may  be  treated  as  usual  by  employing  a  lifting  function.  Furthermore,  we 
define  for  any  system  parameter  vector  p  =  (pi, ...,  pj)  £  V  =  jj[=1  2?,;  a  bilinear  form  o(-,  •;  p)  :  X  x  X  — »  R. 
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and  a  linear  functional  /(•;  fi)  :  X  — >  R  as 


/ 

a(w,  v;  fi)  :=  CLi(w,  v;  /i*),  with 
2=1 


i(w,v;ni)  :=  a((T™ap)  o  7)),  {T™ap)  i(v\ni  Vw,veX, 

i 

f(v;v)  with  :=  f((T?nap)~1(v\ni  oJl)-,^)  Vi;  G  X 


2=1 


(6) 

(7) 


Note  that  the  parameter  dependent  part  of  the  entries  and  the  determinant  of  the  Jacobian  D%  of  the 
mapping  71  are  contained  in  /i,;,  if  necessary.  If  the  material  parameters  do  not  depend  on  the  spatial 
orientation  of  the  component,  the  application  of  l~[ot  to  the  parameter  dependent  vector-valued  functions 
results  in  a  cancellation  of  (7”ro*)^1  in  (6)  thanks  to  D%  =  77°*  G-i  ■  This  therefore  removes  parameters 
which  account  for  the  spatial  orientation  of  the  component  from  the  bilinear  form.  Note  that  for  scalar¬ 
valued  functions  we  obtain  this  cancellation  thanks  to  the  fact  that  T[ot  is  an  orthogonal  mapping.  We  also 
define  an  inner  product  and  the  induced  norm  as 


(v,w)x  :=  J^((V 

2=1 


■map\  —  1 


)  Hn,  °  77),  (7^nap)~  (wjni  °  Ti);m)x 


|x  :=  y/(v,v)x- 


(8) 


Although  both  the  inner  product  (•,  -)a'  and  the  induced  norm  ||  •  ||x  depend  on  geometric  parameters,  we  omit 
the  dependency  on  the  parameter  //  in  the  definition  to  increase  readability.  We  require  that  the  bilinear  form 
a(-,  •;  /i)  is  coercive  and  continuous  on  X  with  respect  to  the  norm  ||  •  ||x  and  define  the  respective  coercivity 
and  continuity  constants  as  follows 


a(v,v,fj,)  a(v,w;n ) 

a[fi)  :=  lm  — ^ ^ —  and  :=  sup  sup  77 — 77 — 77 — 77 — 


IX 


uexwex  \\v\\x\\w\\x 


(9) 


Moreover  we  demand  that  /(-;/x)  is  bounded  with  respect  to  ||  •  ||.y  on  X.  Note  that  this  implies  coercivity 
and  continuity  of  a(-,  •;  p,)  and  boundedness  of  /(•;  /),)  with  respect  to  the  norm  |j  •  ||  y-  Finally,  we  define  the 
coercivity  and  the  continuity  constants  of  a(-,-,/x)  associated  with  the  space  X1^  as 


M,  \  ■  f  a(yW;M)  ,  x  a(v,w,  n) 

a  (H)  :=  ml  — ^ — rr^ —  and  is  (fj)  :=  sup  sup 


vex*  \\v\\*x  -  ■  vex"wex«  \\v\\x\\w\\x 

We  may  now  formulate  the  global  variational  problem  which  reads  as  follows 

For  any  find  ue(/x)  G  X  :  a(ue(n),v;  /i)  =  f(v,n),  \/v  G  X. 

Moreover,  we  consider  the  following  global  "truth"  finite  element  approximation 

For  any  n&V,  find  u{jj)  G  X ^  :  a{u{n),  v;  fi)  =  f(v;n),  \/v  G  . 


(10) 


(11) 


(12) 


Note  that  the  port  compatibility  requirement  (5)  and  our  assumptions  on  a(-,-;/r)  and  /(-;/x)  ensure  well- 
posedness  of  (11)  and  (12).  We  omit  here  and  henceforth  a  super  index  Af  for  functions  in  the  possibly 
high-dimensional  FE  space  Xv"  to  improve  readability. 
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2.2.  Static  condensation 

To  formulate  the  static  condensation  procedure  we  decompose  X ^  into  a  space  corresponding  to  the 
component-interior,  which  we  will  call  henceforth  bubble  space,  and  a  skeleton  space  associated  with  ports. 
First,  we  define  a  bubble  space  on  the  archetype  component  as 

Bm  :=  {v  G  XM  :  v\%  =0,  1  <  j  <  P7}, 

and  for  each  instantiated  component  a  mapped  bubble  space 

Bf  :=  span{77nap(w  o  T~ls),  w  G  BU}.  (13) 

Finally,  we  introduce  a  global  bubble  space  PA"  :=  @[=1  PA7  To  define  the  skeleton  space,  we  choose  the 
coupling  modes  or  interface  functions  ipj  k  G  AA"  as  harmonic  extensions  of  the  port  modes,  i.e. 


(VbA  i  x  —  ^  Vuj  G  B  , 


Xj,k, 

0, 


on  t j, 

on  t y  for  j'  £  j,  ’ 


(14) 


for  1  <  k  <  A/p,  1  <  j  <  P7.  On  instantiated  components  i,  1  <  *  <  /  we  may  then  define  ipi,j,k  = 
‘T^nav{'(pjtk  °  7)_1).  The  global  functions 


{1pi‘  ,j‘  ,ki 
Ip i,j,k ; 
0, 


in  $V, 
in  fli, 

in  \  (fip  U  flj), 


(15) 


form  the  global  skeleton  space 

Snv  :=  i^P.k,  1  <  /c  <  A/p,  1  <  p  <  Pr}- 


(16) 


Next,  we  introduce  for  each  instantiated  component  i,  1  <  i  <  I  &  source  bubble  b{ (yt)  G  B1^ ,  defined  as 
the  solution  of 

a{b{(m),w;in)  =  \/w  £  BM ,  (17) 

and  set  bj(y.i)  =  T™ap  (b{  (y-i)  o7)-1).  Note  that  the  source  bubble  b{(yi)  depends  on  yn,  %  -  1 . ...,/  rather 
than  the  parameter  y  on  the  archetype  component.  Additionally  we  define  bubbles  bij^ifi)  G  B ^  associated 
with  the  coupling  modes  ipj,k  as  the  solutions  of 


a(bi,j,k(^i)  +  tpj,k,w;  yt)  =0,  VroG  B* . 


(18) 


Setting  4>yjlk{^i)  =  h,j,k{Vi)  +  i>j,k  and  )  =  7 7nap(4>i,j,k(ljLi)  °  %  *)>  we  may  introduce  the  global 

functions  d>pj.(y)  G  AA"  as 


Tp,fe(Ai) 


,jf ,k(fJji') )  in  Up, 

in  flj,  ^  7Tp  {(b7)?(/  5J  )}■ 

0,  in  f2  \  (Op  U  ilj), 


(19) 


We  highlight  that  we  have  =  SjyT  ©  B ^ .  The  key  result  is  that  thanks  to  the  definition  of  the  coupling 
modes  (14)  the  decomposition 

XN  =  Bu  ©  span{i/ij  >fe,  j  =  1, ...,  P7,  k  =  1, ...,  TVp},  /  =  !,...,/  (20) 
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is  Jf-orthogonal.  Note  that  we  may  also  employ  a  parameter  independent  bilinear  form  a(-,  ■):  X  x  X  —>■  R 
in  (14)  as  in  domain  decomposition  methods  [25]  or  the  CMS  approach  [4,  5[.  However,  we  emphasize  that 
unlike  in  standard  domain  decomposition  methods  and  the  CMS  approach,  where  the  shape  of  the  component 
does  not  vary,  the  decomposition  on  the  global  system  X1^  =  Sj^v  ©  ZFV"  is  in  general  not  orthogonal  due 
to  the  parameter  dependency  of  the  geometric  transformation  7).  Note  to  this  end  that  in  contrast  to  the 
bilinear  form  in  (18)  the  inner  product  in  (14)  does  not  depend  on  geometric  parameters. 

Finally,  we  introduce  for  each  instantiated  component  a  global  function 

ELi  +  Ef=i  Ekt i  Un iW.kitifajAVi)’ 

0, 

with  unknown  coefficient  Pn.jjyfc^)  G  R.  We  write 

/  I  Pr  A/V 

u O)  =  Ui(m)  =  XX  (w)  +  X  X!  Up,k(lJ.)®p,k  (/*)  (22) 

i—  1  i— 1  p=  1  fc= 1 

and  consider 

For  any  fx  GV,  find  u(fx)  G  X 'X  :  a(u(n),v,  fx)  =  f(v;fx),  \/v  G  Sj^v.  (23) 

Note  that  thanks  to  (17)  and  (18)  problem  (23)  is  equivalent  to  (12). 


fti, 

fl\  Cli 


(21) 


3.  Model  reduction 

As  the  static  condensation  procedure  described  in  §2  is  computationally  very  expensive  we  apply  model 
reduction  techniques.  Within  the  component  interior  we  construct  as  in  the  SCRBE  method  [14]  a  re¬ 
duced  basis  (RB)  approximation  [7]  for  each  interface  and  source  bubble  (§3.1).  Additionally  we  construct 
a  low-dimensional  port  space  by  a  pairwise  training  algorithm  [15]  and  a  subsequent  proper  orthogonal  de¬ 
composition  (POD),  which  is  then  used  to  define  a  port  reduced  SCRBE  (PR-SCRBE)  approximation  (§3.2). 
We  prove  that  the  PR-SCRBE  and  hence  the  SCRBE  approximate  problem  are  well-posed  and  that  the 
approximation  is  stable,  which  is  one  new  contribution  of  this  paper. 

3.1.  Static  condensation  reduced  basis  element  method  (SCRBE) 

For  each  instantiated  component  we  define  the  RB  approximations  b{’N (/n*)  G  Bjf  C  B ^  and  bX  G 

B^k  C  as  the  solutions  of 

a(b{’N (Vi),w,  Vi)  =  Vw  e  Bf  ,  (24) 

bibX^k{ni)  +  4>jik,  w;  m)  =  0,  Mw  G  Bfk,  (25) 

for  1  <  i  <  I  and  1  <  j  <  P7,  1  <  k  <  Af-p.  Hence  each  RB  approximation  b{’N(/Xi )  and  bX  k(/Xi)  is 
associated  with  a  different  RB  space  Bjf  or  B^k  of  dimension  <  N  tailored  to  the  respective  bubble,  where 
N  =  1, ...,  Nmax.  Constructing  these  RB  spaces  with  a  Greedy  algorithm  [26]  ensures  a  rapid  convergence 
of  the  RB  approximation  [8] .  For  further  details  on  RB  methods  we  refer  to  [7] .  We  also  define  b{'N (/q)  = 

T™ap (b{'N (m)  °  X)1)  and 

bJ’N(n)  :=  ^5f;Ar(/Xi), 

2—1 
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(26) 


where  we  extend  b{’N (/q)  by  zero  to  0^  for  i'  ^  i.  Setting  and  <^ifc(Mi)  = 

7jnap k{ni)  °7”_1),  we  may  introduce  the  global  functions  ^^.(m)  €  as 


,*(/*»),  in  fV, 

(/**)>  in  fij,  ,  7TP  =  (*',/)}  (27) 

0,  in  fi  \  (fij/  U  fij), 

and  define  the  SCRBE  approximation  space  Iw(/j)  as 

Xn(/j,)  :=  span{$^fe(/r)  :  1  <  p  <  Pr,  1  <  fc  <  A/p}.  (28) 

We  may  then  consider  the  problem 

Find  Ug(v)  G  :  a(ua(n),v;  n)  =  f{v\n)  -  a(bI[N  n),  \/veXN(/j.),  (29) 

and  define  an  SCRBE  approximation  uG  {^l)  as 

u%(n)=u%{n)+bf’N([i).  (30) 

Here  the  index  G  indicates  that  we  consider  a  (standard)  Galerkin  approximation,  i.e.  that  we  employ  XN  (/j.) 
both  as  a  test  and  a  trial  space. 

Alternatively  we  may  perform  a  Petrov-Galerkin  approximation  and  consider 

Find  UpG{n)  G  XN (^)  :  ci(upG(n),v;  =  f{v\n)  -  a(6/;iv(p),  u; /r),  Vu  €  Sjsv.  (31) 

Then  we  define  the  corresponding  SCRBE  approximation  UpG(fi)  as 

upg( m)  =  +  bf,N  ((i).  (32) 

We  emphasize  that  due  to  the  RB  approximation  the  problems  (29)  and  (31)  are  not  equivalent. 

3.2.  Port  reduction 


Next  we  consider  reduced  port  spaces  [15,  18]  of  dimension  n  <  A/p.  Note  that  in  general  the  dimension 
of  the  reduced  port  spaces  may  vary  within  the  global  system  even  for  the  same  port  type.  However, 
for  the  sake  of  simplicity  and  readability  we  assume  that  all  reduced  port  spaces  within  the  system  have 
the  same  dimension  n.  To  obtain  a  fast  convergence  of  the  port  reduced  approximation  and  hence  an 
accurate  approximation  also  for  n  <C  A/p  we  employ  empirical  port  modes,  constructed  with  model  reduction 
techniques  related  to  RB  methods  [15,  18].  The  key  idea  of  the  empirical  port  mode  strategy  proposed  in  [15] 
is  to  include  knowledge  on  the  solution  manifold  on  the  port  into  the  basis  construction  process  to  obtain  a 
fast  converging  port  reduction  approximation.  It  is  exploited  that  the  behavior  of  the  solution  on  an  interior 
port  Tp  is  fully  determined  by  the  geometric  form,  the  considered  PDE,  the  parameter  values  and  the  behavior 
of  the  solution  on  the  non-shared  ports  of  the  two  components  which  connect  at  Tp.  The  purpose  of  the 
pairwise  training  algorithm  proposed  in  [15]  for  the  construction  of  empirical  port  modes  is  thus  to  explore 
the  solution  manifold  induced  by  the  parameter  dependency  regarding  geometric  and  material  parameters 
and  an  introduced  parametrization  of  the  unknown  solution  behavior  on  non-shared  ports  in  a  systematic 
fashion  to  prepare  the  port  modes  for  all  possible  component  connectivities  and  parameter  configurations 
that  may  be  encountered  in  the  online  stage.  For  further  details  we  refer  to  [15,  18]. 

Now  we  present  the  general  port  reduction  framework  as  introduced  in  [15]  and  prove  well-posedness  of  the 
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PR-SCRBE  approximate  problem. 

First,  we  introduce  a  port-reduced  skeleton  space 

Sn  :=  {'t'p.fc,  1  <  k  <  n,  1  <  p  <  Pr}  C  Sj^fv ,  with  n  <  Afp.  (33) 

Furthermore,  we  introduce  a  port-reduced  SCRBE  approximation  space 

X (p)  :=  span{$^fc(^)  :  1  <  fc  <  n,  1  <  p  <  Pr}5  with  n  <  Af-p,  (34) 

and  consider  the  problem 

Find  u”G(p)  G  X*{p)  :  a(u”G(p),  v;  p)  =  f(v;  p)  -  a(bf’N(p),  v;  p),  Vu  G  X*{p).  (35) 

We  may  then  define  the  Galerkin-PR-SCRBE  approximation  as 

Pr  n 

+  W^g(^)  =  Un,G,p,k(^P,M’  (36) 

p=  1  k= 1 

where  u%(p)  =  Ep=i  and  Un,G,P,M  G  R  denote  the  coefficients  obtained  by 

solving  (35). 

Lemma  3.1  (Well-posedness  of  the  Galerkin-PR-SCRBE  problem  (35)).  The  Galerkin-PR-SCRBE  problem 
(35)  is  well-posed.  Moreover,  for  all  p  G  T>  and  f(p)  G  X  there  holds  the  stability  estimate 

iKoMlIxi^^+^Oll/Wllx,  (37) 

Proof.  Thanks  to  our  assumption  that  a(-,-;p)  is  coercive  on  X  with  respect  to  the  norm  II  -Ik,  the  RB 
approximate  problems  (24)  and  (25)  are  well-posed.  The  mapping  from  a  system  parameter  p  G  T>  to  the 
space  X^ (/x)  is  hence  well  defined.  We  recall  that  functions  in  Xffip)  are  zero  on  ports  where  we  prescribe 
Dirichlet  boundary  conditions  because  Pr  denotes  the  number  of  global  ports  on  which  we  do  not  impose 
Dirichlet  conditions.  The  port  compatibility  requirement  (5)  and  the  fact  that  Sn  C  Sj\fv  then  yield  that 
Xff  (ji)  C  X1^  and  thus  that  the  PR-SCRBE  approximation  setting  is  conformal.  Therefore,  the  bilinear  form 
a(-,  •;  p)  is  coercive  on  Xff  (p)  with  coercivity  constant  (p)  and  continuous  with  continuity  constant  is^(p). 
The  Lax-Milgram  Lemma  hence  implies  well-posedness  of  (35).  Exploiting  equation  (17),  the  coercivity  of 
a(-,-;p)  on  X^ ,  and  the  definitions  of  a(-,-\p)  (6),  f(-;p)  (7),  and  bB-N[p)  (26)  yields 

l|6/;jVk)lk  <  -^rxWfmx'.  (38) 

'>A  (/') 

The  stability  estimate  then  follows  by  testing  with  G{p)  in  (35).  □ 

Again  we  may  consider  the  Petrov-Galerkin  formulation 

Find  u*PG(p)  G  X* (p)  :  a(u*PG(p),v;  p)  =  f(v,  p)  -  a(bf  ’N(p),  v;  p),  Vv  G  Sn.  (39) 

We  may  then  define  the  PR-SCRBE  approximation  as 

pr  n 

un,PG^T)  =  bt'N (A4)  +  ^Pg(P)  =  b*'N{p)  +  E  E  ^nfpG,p,fc(/i)^fc(/i)>  (40) 

P=  1  k— 1 

where  u"PG(p)  =  Ep=i  ELi  H^pg^M^M  and  Un,PG,p,k(T)  e  R  denote  the  coefficients  obtained  by 
solving  (39). 
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Lemma  3.2  (Well-posedness  of  the  Petrov-Galerkin  PR-SCRBE  problem  (39)).  Let  us  assume  that  for  all 
p  £  T>  the  discrete  inf-sup  condition 


inf  sup 

weX ”  (fJ.)  vGSn 


frff  =:  ^  w  > » 

IMMMU 


(41) 


is  fulfilled.  Then  the  Petrov-Galerkin  PR-SCRBE  problem  (39)  is  well-posed.  Moreover,  for  all  p  £  V  and 
f(p)  £  X  there  holds  the  stability  estimate 


un,pc(p) II A'  < 


/TO 


■  X 


1  + 


(p) 


„X 


(p) 


yX 


(p)  J 


\\f(p)\\ 


X' 


(42) 


Proof.  The  fact  that  the  spaces  Xff(p)  and  Sn  have  the  same  (finite)  dimension  together  with  the  Banach- 
Necas-Babuska  Theorem  yield  the  well-posedness  of  (39)  (cf.  [27]).  The  stability  estimate  can  then  be 
obtained  by  using  the  a  priori  estimate  (38)  and  the  discrete  inf-sup  condition.  □ 

Note  that  the  well-posedness  of  (29)  and  (31)  follows  directly  from  Lemma  3.1  and  Lemma  3.2  by  con¬ 
sidering  n  =  M-p . 


4.  Certification  framework 

The  goal  of  this  section  is  to  develop  a  certification  framework  for  the  PR-SCRBE  method  which  provides 
a  rigorous  and  efficient  a  posteriori  error  estimator  for  the  error  j|tt(/u)  —  uff  (/x)||x-  T°  this  end  we  adopt 
techniques  from  the  RB  framework  [7]  to  define  an  a  posteriori  error  estimator  for  the  error  caused  by  the  RB 
approximation  of  the  bubble  functions  in  the  component  interior  (§4.1.1).  To  motivate  our  proposed  concept 

for  an  a  posteriori  error  estimator  for  port  reduction  we  recall  that  the  "truth"  FE  (static  condensation) 

approximation  u(p)  satisfies  a  weak  flux  continuity  across  global  ports;  to  wit  for  two  components  i  and  i' 
which  connect  at  a  global  port  Tp,  np  =  {( i,j ),  there  holds 

fi ('hp./c i  Pi')  u-i ( G  (pi ) .  4/^  ^ ,  pi)  T  f ji  (tkp  ^ ,  pi' )  ap  (up  (pp ) ;  ^ ,  pp )  0,  k  1,...,  J\f*p .  (^3) 

For  the  port  reduced  SCRBE  approximation  we  also  have  weak  flux  continuity  but  only  with  respect  to  the 
reduced  port  space.  To  measure  how  well  the  PR-SCRBE  approximation  behaves  like  the  truth  solution  on 
the  ports  and  thus  to  quantify  the  error  caused  by  port  reduction  we  hence  compute  weak  fluxes  for  the 
PR-SCRBE  approximation  on  all  ports  of  a  component  with  respect  to  the  full  port  space  and  use  the  jump 
of  those  fluxes  across  ports  to  define  an  a  posteriori  error  estimator  (§4.1.2).  To  combine  both  error  estimator 
contributions  we  exploit  the  X-orthogonality  of  the  bubble  and  the  skeleton  space  on  each  archetype  reference 
component  (20).  We  state  the  main  results,  namely  the  a  posteriori  error  estimate  and  the  effectivity  results, 
in  Sections  4.1  and  4.2  and  collect  the  respective  proofs  in  §4.4.  We  employ  approximations  of  the  occurring 
constants  in  the  estimator  and  discuss  in  §4.3  how  such  approximations  may  be  obtained.  We  emphasize 
that  all  results  in  this  section  hold  true  both  for  the  Galerkin  (35)  and  the  Petrov-Galerkin  formulation  (39) 
of  the  PR-SCRBE  method.  Thus  we  omit  the  respective  indices  G  and  PG  in  this  section  and  recall  that 
the  respective  PR-SCRBE  approximations  have  been  defined  in  (36)  and  (40)  as 

Pr  n 

Un  (P)  =  bf’N(p)  +J2Y1  Un,p,k(P)®p,k(P)>  fOT  Un,p,k(p)  G 

p—1  k—1 


li 


4.1.  Formulation  and  analysis  of  a  rigorous  a  posteriori  error  estimator  for  PR-SCRBE  based  on  local  indi¬ 
cators 

4-1.1.  An  a  posteriori  error  estimator  for  the  RB  error 

First  we  consider  error  contributions  due  to  reduced  basis  approximations  of  the  bubbles  defined  in  (24) 
and  (25).  For  each  instantiated  component  i,  1  <  i  <  I  and  1  <  j  <  P7,  1  <  k  <  n  we  define  the  residuals 

r{{w;Pi)  ■=  f{w;  pf)  -  a{b{'N pi),  Vw  G  BM , 
and  fitjtk(w-,Pi)  ■■=  -a(4>j,k  +  bfjk{pi),w-,pi),  Vu>  G  BM . 

Furthermore,  we  introduce  the  corresponding  Riesz  representations  1Z{  (p^  G  B ^  and  Rj.hk(Pi)  G  as 
the  solutions  of 


(1l{ (, pi ),  w)x  =  r{ (w;  pi),  \/w  G  BM , 

(R-ip  ,k(PJi)  ?  X  —  Pi) ,  Ww  G  B 


(44) 


We  may  now  employ  these  Riesz  representations  to  define  an  a  posteriori  error  estimator  A N (p)  for  the 
RB  error  as 


A  N(p) 


with  A  ?(p) 


'-'app'~'2,app 

& app (aO 


(£(A  f(M))2)172 


pT 


ffdlGilDI  detG, 


|-i/2| 


R-i  ( Pi )  +  Bn,Ili(j),k(p)Ri,j,k(Pi)\\ X  ■ 

3  = 1  k—1 


(45) 


We  emphasize  that  the  inner  product  employed  in  (44)  does  not  depend  on  (geometric)  parameters.  Therefore 
the  RB  error  estimator  (45)  is  based  on  information  on  the  error  on  the  (archetype)  reference  domain 
and  not  on  the  mapped  domain  fl;  in  the  global  system.  Note  that  this  is  compliant  with  the  standard 
RB  methodology  and  error  estimator  for  parameter  dependent  domains  O  [7].  However,  the  parameter 
dependency  of  Qi  is  reflected  by  the  factor  g(||Gi||)|  detGi|-1/2,  on  which  we  will  elaborate  below.  Note  also 
that  due  to  the  fact  that  we  employ  a  different  RB  approximation  for  each  bubble  associated  with  a  coupling 
mode  (cf.  §3.1)  we  have  to  consider  a  linear  combination  of  Riesz  representations. 

We  also  propose  a  second  error  estimator  A N,1(p)  for  the  RB-induced  error  which  is  computationally  more 
feasible  than  A N  (p)  especially  in  terms  of  storage  (see  §5).  We  define  A N,1(p)  as 


rX  rX  1 

C„™,C9  .  N1  2\  I/2 


A N’\p)  W) 

°^app\H')  i=1 

P ^  n 

for  A  ^\p)  ~  g(\\Gi\\)\detGi\~1/2  (\\K{  (p^y  +  J2J2\Un,npj),M\\\Ri,3,k{Pi 


(46) 


IX 


3= 1  fc=l 


We  remark  that  one  vital  improvement  to  previous  contributions  [14,  15]  is  that  both  A N (p)  and  AN'1  (p) 
yield  an  estimator  for  the  X-nonn  of  the  error  (see  Proposition  4.2  and  Corollary  4.3)  and  that  their  validity 
does  not  rely  on  the  quality  of  the  RB  approximation.  We  also  highlight  that  both  A N(p)  and  AN'1(p) 
consist  of  local  indicators  allowing  for  an  adaptive,  component-wise  choice  of  the  bubble  spaces. 

Regarding  the  constants  in  (45)  and  (46),  we  note  that  aapp(p)  is  an  approximation  of  a^(p)  or  a(p)  and 
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we  set  cx  :=  sup„eX.v  |M|ffi(n)/IMU  and  cf  :=  sup„eXJv  |MU/|M|ffi(n)  and  denote  with  c*p  and  c£opp 
the  respective  approximations.  We  remark  that  we  could  alternatively  define  an  a  posteriori  error  estimator 
for  the  RB  error  by  omitting  both  cxpp  and  cx app.  However  adding  those  constants  improved  the  overall 
performance  of  the  PR-SCRBE  estimator  significantly  in  our  numerical  experiments.  Finally,  we  remark 
that  the  definition  of  the  function  g  :  R+  — »  R+  in  <7(||Gj||)  depends  on  the  chosen  inner  product 
For  the  H[  {il) rfr -sem i-norin  we  have  for  instance  g(||Gj||)  =  ||Gj||  and  for  the  full  W  (0)dr-norm  we  obtain 
g{ 1 1 G* 1 1 )  =  1 1 Gj||  +  1.  For  other  choices  of  the  inner  product  (-,-)a  the  precise  definition  of  g(||Gj||)  follows 
from  straightforward  calculations  exploiting  the  transformation  theorem. 


^.1.2.  An  a  posteriori  error  estimator  for  port-reduced  approximations  based  on  conservative  fluxes 

Next,  we  address  the  error  caused  by  port  reduction.  To  this  end  we  analyze  the  jump  in  the  flux  of 
the  PR-SCRBE  approximation  ux  (p)  over  the  ports.  As  motivated  at  the  beginning  of  this  section  we 
wish  to  consider  a  weak  flux  and  hence  propose  to  use  a  conservative  flux  defined  according  to  Hughes  et 
al.  [19]  as  a  flux  measure.  The  main  result  in  [19]  is  that  by  employing  a  conservative  flux  the  continuous 
Galerkin  method  is  locally  conservative  with  respect  to  subdomains  consisting  of  a  union  of  grid  elements. 
The  work  is  based  on  earlier  results  (see  for  instance  [28,  29])  which  obtain  global  conservation  in  spite  of  the 
presence  of  Dirichlet  boundary  conditions  by  post-processing.  Note  that  the  conservative  flux  is  a  variational 
approximation  of  the  exact  flux  over  the  boundary  of  the  component.  We  adopt  the  concept  proposed  in  [19] 
to  the  PR-SCRBE  setting  by  defining  a  conservative  flux  for  1  <  *  <  I  as  the  solution  of 


=  fi(v;fJ-)  ~  ai(ux{p),v;p),  VueSVp,  (47) 

where  E i  =  Ui<j<p-r  7 i,j-  Note  that  thanks  to  our  mutual  disjoint  port  assumption  and  the  definition  of  the 
coupling  modes  no  post-processing,  i.e.  computation  of  a  global  conservative  flux,  due  to  Dirichlet  boundary 
conditions  as  described  in  [19]  is  required.  Note  also  that  provided  =  1  e  Sj\fv  we  may  verify  mass 
conservation  and  that  Hfffp)  is  indeed  the  conserved  total  flux  along  E,;.  Due  to  the  mutual  disjoint  port 
assumption  problem  (47)  decouples  over  different  ports  and  we  may  compute  the  conservative  flux  separately 
for  each  port.  Thus  we  have  Hx^{p)\li:j  =  where  the  latter  is  defined  as  the  solution  of 


fi  (ifi.j.k'-!  P)  ai(un  (lO)  A^)j 


1  <  j  <  P1, 1  <  k  <  Mv. 


(48) 


We  may  then  define  the  jump  of  the  conservative  flux  across  a  global  port  Tp  as 

[Hn]P(x-,p)  :=  Hxid(x-,  p)  +  p),  where  ttp  =  {(i,j),  (*',/)}• 


(49) 


If  a  local  port  7 ij  lies  on  the  part  of  1  where  Neumann  or  Robin  boundary  conditions  are  prescribed  we  set 
[. Hx]p(x ;  p)  :=  Hfft  Ax ;  p).  Based  on  that  we  introduce  an  a  posteriori  error  estimator  for  port  reduction 


An(/i) 


d'  cx 

^app^app 

^appi'fd) 


Pr 


®AP») 


2D/2 


for  A P(p)  :=  \\[Hx]p(p)\\L2irp)dr. 


(50) 


Here  :=  maxc(  app  and  c\  app  is  an  approximation  of  c\  :=  sup„eX^  ||v|||2(I].)/||^||^i(n.)  for  component 

i.  We  also  introduce  ct  :=  rnaxc*. 

iei  1 

We  highlight  that  in  the  spirit  of  standard  residual  based  a  posteriori  error  estimators  for  the  FEM  A n(p) 
consists  of  local  indicators  A([(//).  Below  we  show  a  local  effectivity  result  (see  Proposition  4.5)  for  the  local 
indicators.  Hence  they  may  be  used  within  an  adaptive  PR-SCRBE  scheme  for  an  efficient  choice  of  local 
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port  spaces  on  each  global  port.  Note  that  this  is  one  new  contribution  of  our  proposed  error  estimator,  as  the 
error  estimator  proposed  in  [15,  18]  is  based  on  a  (global)  non-conforming  version  of  the  Schur  complement 
matrix. 

Furthermore,  we  emphasize  that  we  only  compute  the  flux  on  the  ports  of  each  component  and  hence  avoid 
the  computation  of  a  residual  on  the  full  global  port  space  as  in  the  CMS  approach  [5].  Moreover,  thanks 
to  our  assumption  that  7),  i  =  1, ...,/  when  applied  to  a  port  is  a  rigid  body  motion,  the  port  modes  are 
L^yi^-orthonormal,  i  =  1 ,...,/,  j  =  1  ,...,P7,  and  the  computation  of  the  conservative  flux  reduces  in 
this  case  to  the  assembling  of  the  residual.  If  we  allow  geometric  mappings  7),  i  =  1 that  deform  the 
port,  the  port  modes  are  in  general  not  L2(7,j)-orthonormal  anymore,  and  the  computational  costs  for  the 
computation  of  An(p)  are  therefore  higher.  The  effect  on  the  computational  costs  are  further  discussed  in 
Subsection  5.2.  Note  that  by  restricting  the  test  space  in  (47)  to  a  subset  of  v  a  hierarchical,  non-rigorous 
estimator  may  be  computed  at  lower  cost.  As  in  many  cases  the  coefficients  of  the  truth  static  condensation 
FE  approximation  for  the  high  modes  are  very  small,  we  expect  that  also  this  cheaper,  non-rigorous  estimator 
would  provide  a  very  good  estimate  of  the  error  induced  by  port  reduction. 


^.1.3.  Properties  of  the  decomposition  X ^  =  B1^  ©  Sjy 

Recall  that  the  RB  error  estimator  (45)  is  based  on  the  dual  norm  of  weighted  residuals  with  respect  to  the 
bubble  space  B1^  and  that  the  PR  error  estimator  (50)  rests  upon  the  jump  of  the  conservative  flux  using  the 
skeleton  space  S/s  as  a  test  space.  To  derive  an  a  posteriori  error  estimator  for  the  error  u(p)  —  uff  (p)  £  X1^ 
in  the  X-norm  we  thus  require  estimates  of  the  type  ||'P|jjc  <  Ci||u||x  and  ||6||x  <  C2||v||x  for  all  4/  £  Sj^p, 
b  £  B*  and  4'  +  b  =  v  £  A‘v" .  The  two  essential  ingredients  to  derive  such  estimates  are  the  decomposition 
XA"  =  ©  Sj\f  and  the  fact  that  the  decomposition 

XN  =  Bm  ©  span{4fe,  j  =  l,...,P\k=l,...,Uv},  i  =  l,...,I  (51) 

is  X-orthogonal  (see  §2.2  for  further  details).  We  collect  the  results  in  the  following  proposition. 

Proposition  4.1  (Stability  of  the  decomposition  X-v"  =  B1^  ©  Sjvp).  Under  the  assumptions  of  Section  2.1.2 
on  the  mapping  71  the  decomposition  X^  —  B ^  ©  5V  stable  in  the  sense  that  for  all  v  £  X1^ ,  expressed 
as  v  =  4'  +  b  with  4*  £  5v  and  b  £  X ^  we  have 

mh+Mxi  <27(l|G'i||)2ff(l|Gl-1||)2|HI^,  i  =  l,  (52) 

for  constants  0  <  ff(||Gi||),  (?(||G,“1||)  <  00.  Moreover,  the  following  estimates  hold  true 

||^lk<fl(||G<||)5(||G71||)|H|A1  and  \\b\\Xi  <  ff(||Gi||)fl(||GT1||)||t,||A1,  i  =  l,...,L  (53) 

If  the  transformation  71  admits 

(ipi,j,ki  b)Xi  =0  Vb£Bf,  i  =  1, ..,  7,  (54) 

we  have  the  improved  result 

Wnh  +  lMx^M^  i  =  (55) 

Note  that  for  instance  pure  rigid  body  transformations  allow  for  (54). 
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4-1-4-  An  a  posteriori  error  estimator  for  PR-SCRBE 

Based  on  the  results  above  we  may  now  introduce  an  a  posteriori  error  estimator  A(yrt)  for  the  error  caused 
by  the  PR-SCRBE  approximation  as 

A  (p)  =  An{p)  +  AN(p).  (56) 

The  estimator  A(p)  is  a  rigorous  and  robust  bound  as  stated  in 

Proposition  4.2.  Let  the  approximations  of  the  constants  fulfill  cx / 'a(p)  <  cxpp/aapp(p),  cx  (p)  <  cxapp(p) 
and  c*  <  c*pp.  Then  the  error  estimator  A(p)  =  An(p)  +  AN(p)  with  An(p)  and  AN (p)  defined  in  (50)  and 
(45)  satisfies 

\\u(p)-u%(p)\\x<A(p).  (57) 

Based  on  the  second  RB  error  estimator  AN,1(p)  we  define  a  second  a  posteriori  error  estimator  for  the 
error  caused  by  the  PR-SCRBE  approximation  as 

A\p)  =  An(p)  +  AN’\p)  (58) 

for  which  we  obtain  the  following  result. 

Corollary  4.3.  Let  the  approximations  of  the  constants  fulfill  cx/a(p)  <  cxpp/aapp{p),  cx(p)  <  cxapp(p) 
and  <*  <  c^pp-  Then  the  error  estimator  A1(p)  =  An(p)  +  AN,1(p)  with  An(p)  as  in  (50)  and  AAr,1(/x)  as 
defined  in  (46)  satisfies 

\\u(p)-ux(p)\\x<A1(p).  (59) 

Note  that  the  structure  of  A(p)  and  A1(/r)  allows  us  to  balance  the  RB  and  the  PR  error  contributions. 
Moreover,  if  the  RB  error  is  small,  we  may  employ  only  An(p)  as  a  (non-rigorous)  a  posteriori  error  estimator, 
whose  computation  reduces  in  general  to  the  assembling  of  the  residual  in  (48). 


4-2.  Analysis  of  the  effectivity 

In  this  subsection  we  derive  bounds  for  the  effectivities  of  the  local  indicators  A(v(/.i)  (45)  and  A^(p) 
(50).  Based  on  that  we  conclude  that  also  the  effectivities  of  AN (p),  An(p),  and  thus  A(/i)  are  bounded. 
First,  we  study  the  effectivity  of  AN (//,),  i.e.  we  consider 


pN(p)  :=AN(p)/\\u(p)-u"(p)\\x, 


and  obtain  the  following  result. 

Proposition  4.4  (Effectivity  of  AN(p)).  The  local  error  indicators  Ax (p)  satisfy 

Af(/u)<ZA(/r)5(||Gi||)5(||G-1||)||u(M)-^(/z)||xj. 


As  a  result  we  obtain  the  following  bound  for  the  effectivity  rjN  (p) 

'"(p) 


VN{t)  < - f^capPc2,apP  max  o(||Gi||)5((||G.  1\\). 

(Xappxf-l)  z— 


(60) 


(61) 


(62) 


We  recall  that  we  may  omit  cxppcxapp.  Thus  if  77,  *  =  ■■■,!,  is  only  a  composition  of  a  translation  and 

a  rotation  and  hence  ||Gj||  =  1 1  G~ 1 1 1  =  1  we  recover  the  standard  effectivity  result  for  the  RB  error  estimator 
based  on  the  dual  norm  of  the  residual  [7].  Note  also  that  against  the  background  of  Proposition  4.1  the 
constant  depending  on  Gi  in  (62)  reflects  how  much  the  geometrical  transformation  affects  the  decomposition 
©  SVP  in  terms  of  a  deviation  from  orthogonality. 
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We  remark  that  we  cannot  obtain  an  effectivity  result  for  AN’1(pf 


As  a  second  step  we  derive  a  bound  for  the  effectivity  of  A n(p) 


Vn(p)  ■=  An(/z)/|| u{p)  -  (m)|U, 


(63) 


as  demonstrated  in  the  following  proposition. 

Proposition  4.5.  We  obtain  for  the  local  error  indicators  A ^(p)  that 

A £(/x)  <  ^(/r)ch/i“1/2C,p(/u)||u(/u)  -  (/u)lki®xi,,  for  ttp  =  {( i,j ),  (*',/)},  (64) 


wftere  |M|xi0x.,  :=  (IMIx,  +  IMIiv)1/2>  <??(/*)  ==  1)1  det  G-i | 1/2  +  ct-yff(||Gi(1||)|  det  GV|1/2,  c^- 

is  the  constant  in  the  trace  theorem,  h  denotes  the  mesh  size  of  the  port  mesh,  and  Ch  is  the  constant  in 


the  inverse  estimate 
satisfies 


|ffi/2(rp)dr  <  chh  1/2||u||L2(rp)d.  [30]  for  v  G  Moreover,  the  effectivity  rjn{p) 


Vn(p)  < 


( 

O^avp  (aO 


oC-app'/P^Chh  1/2  max  Cp(p). 


P=i,...,Pr 


(65) 


We  emphasize  that  in  actual  practice  the  dependence  of  r]n{p)  on  h  is  weak,  as  will  be  demonstrated  in 
the  numerical  experiments  in  §6.  Alternatively  we  may  consider  J/1  / 2 -ort lionorm al  port  modes,  obtained  by  a 
lifting  procedure  and  a  corresponding  conservative  flux  expressed  as  a  linear  combination  of  these  modes.  In 
this  way  we  can  derive  a  bound  for  the  effectivity  r]n(p)  which  is  independent  of  the  underlying  finite  element 
discretization.  For  further  details  on  the  definition  and  computation  of  H1  /2-orthonormal  port  modes  we 
refer  to  Appendix  A. 

Finally  we  introduce  the  effectivities  of  A(p)  and  A  1(p) 

r]{p)  :=  A(p,)/\\u(ii) -u*{p,)\\x  and  ^(p)  :=  A1(p)/\\u(p)  —  u%{p)\\x  (66) 


and  remark: 

Corollary  4.6  (Effectivity  of  A(/x)).  The  effectivity  r/(p)  can  be  bounded  as  follows: 


n(p)  < 


^(p) 

aapp(p)l'apPV'2’a™i=^,i 


X  I  X 


max  o(||Gi||)5((||G.  1 1|)  +  c*  s/P^chh  1/2  max  Cp(p))  .  (67) 

2=1,. ..,1  p=l,...,Pr  J 


Proof.  The  bound  (67)  is  a  direct  consequence  of  the  definition  of  A(/i)  (56),  Proposition  4.4  and  Proposition 
4.5.  □ 


].3.  Estimation  of  constants 

This  subsection  is  devoted  to  the  discussion  of  various  possibilities  to  determine  approximating  constants 
aapp(p),  c^pp,  c]fpp  and  c*app  of  the  constants  in  the  a  posteriori  error  estimators  A(p)  and  Ax(p). 

In  some  cases  a  lower  bound  for  a(p)  and  thus  cX  (p)  can  be  derived  analytically  as  for  instance  for  heat 
conduction  problems.  In  general  we  propose  to  use  ct^(p)  =  inivGX^(^)a{v^vi  p)/\\v\\x  which  may  be 
determined  by  computing  for  a  given  parameter  p  G  T>  the  smallest  eigenvalue  of  the  eigenvalue  problem: 
Find  the  set  of  eigenvalues  and  eigenvectors  (A^ ( p),v % (p)),  where  (p)  G  K+  and  (p)  G  X]f  (p)  satisfy 

a(Vn{p),w-,p)  =  \%(p){v%(p),w;p)x  Vru  G  (p).  (68) 
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As  the  left  hand  side  in  (68)  yields  the  Schur  complement  matrix,  this  requires  the  assembling  of  the  inner 
product  on  the  right  hand  side  and  the  solution  of  the  eigenvalue  problem  (68).  Alternatively  we  may  compute 
the  smallest  eigenvalue  of  the  Schur  complement  matrix  which  we  expect  to  be  a  good  approximation  (but 
not  necessarily  a  lower  bound)  of  qA  (yrt)  thanks  to  the  minimax  principle. 

It  remains  to  justify  the  choice  ctapp(fi)  =  qA  ( /-/) •  Note  that  the  derivation  of  a  rigorous  lower  bound  of 
gA  (/z)  is  beyond  the  scope  of  this  paper.  However,  we  expect  a ^  (/z)  to  be  a  good  approximation  of  oA”(/z) 
for  the  following  reasons.  We  denote  by  A  and  M  the  matrices  associated  with  the  left  and  right  hand 
sides  of  the  "truth"  eigenvalue  problem:  For  given  /z  G  V  find  (A -^(/z),  i>(/z))>  such  that  A^(/i)  G  R+  and 
v(n)  €  X1^  satisfy  a(v(n),w;  /z)  =  A w)x  for  all  w  G  X1^ .  Moreover  the  subindices  P  and  i  refer  to 
the  degrees  of  freedom  associated  with  the  ports  in  the  system  and  the  interior  of  component  i,  respectively. 
Exploiting  the  decomposition  X ^  =  S yvy,  and  suppressing  the  dependency  of  /z,  the  Schur  complement 

system  for  the  “truth”  eigenvalue  problem  for  the  smallest  eigenvalue  (n)  =  qA(/z)  then  reads  [31]:  Find 
(A^"(/z),  up(/z))  G  M+  x  Mp  such  that 

S(\Y)up(fi)  =  0,  where 

(69) 

S( Xff)  =  ( APP  -  X ?MPP)  -  [(aIp  ~  A Ym^KAu  -  X«  Mn)-\AiP  -  X* MiP)\  , 

i= 1 


where  u p(/n)  denotes  the  coefficients  of  the  P  =  PTJ\f-p  degrees  of  freedom  on  the  ports  in  the  whole  system. 
We  observe  that  (Au  —  A J^Mii)~1(AiP  —  Xff  MiP)  corresponds  to  the  component-local  problems:  For  given 
{m,cr)  G  Vi  x  [0 ,aXr(ni))  find  S?Afc(/Zi)  G  BM  such  that 


a(bZj,k(tJ'i)>w'i  fa)  -  G  BN ,  (70) 

for  i  =  1 j  =  1  ,...,P7,  k  =  1,...,A fp.  Here  off  (/z*)  :=  inf a(D,  ti; /Zj)/||'0||^.  .  and  o  parametrizes 

the  smallest  system-level  eigenvalue  Xff  (/z)  =  qA”(/z).  Note  that  a  similar  approach  has  been  proposed  in 
[32],  We  emphasize  that  we  request  o  G  [0, off (m))  and  hence  qA(/z)  <  off  (/Zj)  to  obtain  well-posedness  of 
(70).  Otherwise  we  define  aapp(lP)  '■=  minj=i;.../  offil^i)- 
Simple  calculations  yield  the  estimate 


-  hj,k{Vi),bXk(m)  ~  ffjAlM);  Mi))1/2  <  v{aff{ni))  1/2  !v^l+Cr  11^. 


di  (Mi)  - 


’jAIIxj/I’ 


for  i  =  1, j  =  1,...,P7,  k  =  1,  A/p,  Cff(lM)  ■=  sup^^  sup^^  a(l),ui;  Mi)/||£||;f.  Jw|lA;ia>  and 


frjj,fc(Mi)  as  defined  in  (18).  Hence  for  qA (m)  <C  aff(ni)  we  can  approximate  S( X^)  in  (69)  with  Sapp( Xi  app)  = 
(■ APP  -  X ffppMpp)  -  E[=1  (AjP  -  XffppAffp)AffAiP,  which  justifies  the  choice  aapp( /z)  =  (m) in  the  case 
of  no  RB  error  and  n  =  A/p.  We  emphasize  that  for  problems  where  qA"(m)  depends  on  the  geometry  of  O 
as  in  the  case  of  linear  elasticity  we  expect  a-A^(/z)  off  (/z*),  which  is  verified  in  the  numerical  experiments 
in  Section  6. 

Thanks  to  the  rapid  convergence  of  the  RB  approximation  [8]  we  expect  also  that  6  A  fe(Mi)  is  a  good  approxi¬ 
mation  of  fe(/zz).  Moreover,  due  to  the  expected  good  approximation  behavior  of  the  empirical  port  modes, 
demonstrated  in  Section  6  and  [15,  18],  we  expect  a  fast  convergence  of  qA (/z)  to  the  smallest  eigenvalue  of 
Sapp{XPapp)  [31,  33],  which  finally  justifies  the  usage  of  cff  ( /z )  as  an  approximation  of  ot/^(n). 


As  we  apply  the  trace  theorem  locally  for  each  instantiated  component  in  the  proof  of  Proposition  4.2 
we  may  exploit  the  well-known  result  that  the  optimal  constant  c\,  i  =  1, is  the  largest  eigenvalue  of  the 
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generalized  eigenvalue  problem:  Find  the  set  of  eigenvalues  and  eigenvectors  (A (/x),u(/t)),  where  A(yrt)  £  R+ 
and  v(n)  €  satisfy 

(v{n),w)L 2(E;)  =  {07nap)~1  °  %),  {Tfnap)~1{w  o  7i);  Hi)L2t£\ 

,  (71) 

=  Kn){(VWP)  1Wv)oTi),(J™ap)  =  \(^){v(n),w)Hi{ni)  VweXf. 

Note  that  due  to  the  parameter  dependency  of  the  norm  standard  techniques  as  the  successive  constraint 
method  [7]  are  in  general  not  applicable.  In  some  cases  as  for  simple  dilations  we  can  determine  analytically 
for  which  parameter  we  obtain  an  upper  bound.  For  general  geometric  transformations  we  propose  to  solve 
the  eigenvalue  problem  for  the  parameter  independent  norms  on  the  archetype  reference  component,  i.e.  we 
consider:  Find  the  set  of  eigenvalues  and  eigenvectors  (A,  v),  where  A  €  M+  and  v  £  X^  satisfy 

{v,w)L-*(±)  =  A(D, «>)#!(„)  Vw£XX.  (72) 

Then  we  proceed  as  in  Proposition  4.1  to  arrive  at  the  following  estimate  for  the  trace  constant 

cJ<ct(||Gj||  +  l)|detGj|-1/2=:<aRp.  (73) 

where  cl  is  the  largest  eigenvalue  of  (72).  Note,  that  it  is  essential  to  compute  a  rather  sharp  bound  for  c\, 
i  =  1, as  this  constant  balances  the  contributions  of  the  RB  and  PR  error  estimators. 

Finally  we  address  the  approximations  cAp  and  cA app,  where  we  focus  on  the  more  involved  case  cApp.  For 
many  PDEs  estimating  cA  requires  the  application  of  Friedrichs  inequality.  To  compute  an  approximation 
of  the  constant  in  this  inequality  we  propose  to  proceed  as  in  [34]  and  decompose  O  in  non-overlapping 
subdomains  on  which  we  can  derive  estimates  for  the  local  constants.  We  use  the  sharp  bound  diarn(0,)/7T 
for  the  constant  in  Poincare’s  inequality  for  bounded,  convex  domains  [35],  solve  the  eigenvalue  problem 
corresponding  to  Friedrichs  inequality  on  the  archetype  reference  component,  and  use  again  Proposition  4.1 
to  arrive  at  a  bound  for  the  constant  in  Friedrichs  inequality  with  respect  to  Q.  For  further  technical  details 
we  refer  to  [34]. 

4-4-  Proofs 

Proof  of  Proposition  4-1-  By  expressing  v  £  X A”  as  v  =  'Tfnap(v  °  7”  1),  with  v  £  W  and  invoking  the 
transformation  theorem  we  obtain  for  all  v  £  X^ 

IMI-y,  <  5(l|G"1||)|detG,:|1/2||u||^  and  ||t)||1-  <  5(||G2;||)|  det  G.i|'1/2||u||Xi-  (74) 

Then,  for  v  £  Wv",  expressed  as  v  =  \F  +  b  with  'F  £  S_\fp  and  b  £  ,  we  exploit  (74)  to  obtain 

ll*lk  <5(l|G-1||)|detGi|1/2||4,||^  and  \\b\\Xi  <  5(I|G-1||)|  det  Gi|1/2||S||*. 

As  4/  £  Sj*f  ,  we  may  write  \F  as  a  linear  combination  of  the  coupling  modes  V’i.fc)  *  =  1, . . . ,  A,  j  =  1,  ...,P7, 
k  =  1, ...,  Afp.  Thanks  to  (14)  also  \F  thus  fulfills  ('F,  b)x  =  0  for  all  b  £  B1^  which  implies 

m%  +  ml  =  m2x-  (75) 

As  a  direct  consequence  we  have 

11*11*  <  INI*  and  WHx  ^  INI*-  (76) 
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Next,  we  apply  again  (74)  to  conclude  that 

||^||xi<5(||Gr1||)|detGi|1/Vl|Gi||)|detGr1/2||^lk,  =ff(||Gr1||)5(||Gi||)||V||xi 

and  \\b\\Xl  <5(l|G-1||M||Gi||)||z;||xI. 

Squaring  and  summing  up  yields  (52). 

Finally,  we  assume  that  7 1  admits  (54).  Expressing  T  as  a  linear  combination  of  the  mapped  coupling  modes 
ipi,j,k,  *  =  1, ...,  /,  j  =  1, ...,  P':  ■  k  =  1,  and  exploiting  (54)  yields  (\F,  b)xi  =  0  and  thus  the  claim.  □ 

Proof  of  Proposition  4.2.  As  ®  6V„  ,  we  may  write  each  v  €  as  v  =  b  +  with  b  €  B x  and 

^  e  Stfp.  Then  we  have 

a(u(p)  -  ux (p),v;  p)  =  f{v;  p)  -  a{ux (p),  v;  p)  =  f(b;  p)  -  a(ux (p),  b;  p)  +  /OF;  m)  -  a(u % (p),^\p).  (77) 

N _ v^_ _ ✓  S _ v*_ _ ✓ 

I:=  II:= 

For  the  first  part  I  we  invoke  the  definition  of  the  R.iesz  representations  (44)  and  the  bilinear  (6)  and  linear 
forms  (7)  to  obtain 


I  P"1  n 

1  =  ~  ai(bt'N(Pi)  +  E  E  c/^ni(j),fc(M)^>fc(Mi),  b;  Pi)) 

i—1  j= 1  k= 1 

I  pi  n 

=  (pi)  +  li(j),k(p)'R-i,j,k(pi)ib) X 

i—1  j= 1  fc— 1 

J  pi  n 

<  £  \\M (Mi)  +  E  E  KniU),k {p)Ki,jM\\x\\Hx- 

i—1  j—1  k—1 

Exploiting  that  functions  in  B-f  vanish  on  global  ports  (see  (13))  the  second  part  II  can  be  estimated  as 
follows: 


H  ~  E  (Hn,i(k)’  ^ )L2C£j)dr  ~  Ed^nkW.  ^ )l2(Tp )dr 


i-1 

Pr 


P=  1 


I  P~< 


<(Y\\{HnW\\lnrpyr) (EE  II* 

P=1  i=l  j=l 

=  (Eii[Hawiih(r,,®1/a(Ef:ii»e.hu) 

p— 1  i=l j— l 


2  \l/2 

i2(7«ddG 


.1/2 


pi  /  pi 

T/2t„ t  ST  II  , ,||2  \l/2 


<  (EllMp(M)l|2p2(rp)^)  7  (^El^ll^fn,^)  7  <  (EIIMp(m)IIWp)E^ 1/Vct||u|U- 

p— 1  7—1  p— 1 

Proposition  4.1  then  leads  to 

pr 

T/2 


a{u{p)-u% {p),v\p)  <  cx  [c*(E  ll[^]p(M)lli2(rp 


P=1 


J  Pi  n 

+  ^(E^limetGr'llE (Mi)  +  E  E  <niO),fe(M)^i,/,fc(Mi)ll^)172]  IMk- 

i=i  fe=i 


i=l 
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Exploiting  the  coercivity  of  o(-,  •;  /_/)  with  respect  to  the  space  yields  the  claim.  □ 

Proof  of  Proposition  4-4-  For  notational  convenience  we  introduce 

P7  n 

Mpd  ■■=  M  (pt)  +  E  e  ^ 

j=i  fc=i 

and  set  7 Zi(m)  =  T™ap(R-i(pi)  o7)_1).  Exploiting  the  definition  of  rt(^i)  in  (21), (22),  the  definition  of  the 
bubbles  in  (17)  and  (18),  and  Proposition  4.1  we  may  estimate  on  each  instantiated  component 

p 7  n 

\\M (Pi)  +  E  E  Un,niU)A»)KijM\\2x 

j= i fe= i 

P7  n 

=  fi(J^-i(Pi)’i  Pi)  —  ai(b{  (Pi)  +  ^  1  ^n.IIj  f  i).k(P)<t>i.i,k  (Pi)i  T^'i(Pi)’i  Pi) 

3= 1  fe=l 

P7  A/'t.  P7  n 

=  ai(^i  (l4/)  "b  ^  ^  ^  ^  Uui(j),k(p)(t)i,j,k(Pi)  —  frj  (/^i)  "b  ^  An.U4<i).k(P)Cl>i.i.k(pi)^'^'i(Pi)'  Pi) 

j— 1  /c— 1  J  =  1  /c— 1 

=  ai(u(p)-Un  (p),Ki(pi)-,m)  <  v" (ii)\\u(ii)  -  u% (riWxMKiin^Wxi 

<  ^(p)\\u(p)  -  (M)lk  «7(||G-1||)|detGi|1/2||^(Mi)llx- 

Squaring  and  summing  up  yields  the  claim.  □ 

Proof  of  Proposition  4-5-  We  consider  the  jump  of  the  conservative  flux  [Hff]p(p)  across  an  arbitrary  global 
port  rp,  7 tp  =  {(i,j),  (i',j')}.  Thanks  to  (49)  we  may  introduce  an  extension  d\p[Hff]p{p)  =  jAk=i  Bk(p)^P,k 
of  [Hff]p(p)  with  coefficients  Qk  €  R.  By  expressing  y{p[H^]p(/j,)  as 

]n iU)(p)  =  rP(\.(j)KW]n,(3)Mo7r1) 

on  O,,  exploiting  the  weak  flux  continuity  of  the  “truth”  FE  solution  u(jj)  across  global  ports  and  using 
Proposition  4.1  we  obtain 

mX(p)\\2mrpy, 

=  Mttp[Hn]p(p);P)  -  aAn  (p),ttp[Hn]p(p);P)  +  fi’  ffip  Wu  ]p(p)i  P)  ~  «*'  («n  (p)  ^[Hu]p(p)\ P) 

=  a,i{u{p)  -  u % {n),m p[hAp(p);p)  +  o-i’  (u(p)  -  u* (p),VKp[H™]p(p)-  p) 

<  ^(p)\Hp)  -  uAp)\\X.exA\\KP[HX(p)\A  +  \mHX(p)\\xP) 

<  (p)\\u(p)  —  un  (p)\\xi®xi,(g(\\Gf1\\)\  detGj|1/2||fHni(j)[ff^r]ni(j)(/i)||^ 

+  5(||GJ7l)|detGd1/2||fHn,,aqMn!,ao(/i)b). 

As  91p[i7,^r]p(/r)  €  SVP  me  may  exploit  (14)  to  obtain  that  94n,(j)[77^]ni(J)(A4)  is  the  harmonic  extension 
of  [Hff]n.^(n)  on  Cl,  where  [H^]ni(j)(p)  =  'T™‘ap  ([Hff]n.^(p)  o7)_1).  Thus  we  may  apply  a  trace  theorem 
[36]  to  obtain  : 

||[Jff9f]J,(Ju)||i2(ri>)^  <  ^(^[[^(m)  -  ^(^)IUi©^,(ct,)^(||G71||)|  detG,|1/2||[^]ni0.)(At)||^l/2(^)^ 

+  cP,i'fl,(l|Gi,1H)|  det  Gp  | 1/2 1|  [Hff]n^ q,) (/r) || jyi/2(^,)<ir ). 
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Thanks  to  our  assumption  that  7)  when  applied  to  a  port  is  a  pure  rigid  body  motion  and  the  port  compat¬ 
ibility  requirement  (5)  we  can  further  estimate  : 

where  Cp(/x)  =  ct'),;c/(||G“1||)|  detG;!1/2  +  ct'!,/gr(||G^1||)|  AetGif1^2.  By  applying  the  inverse  estimate 
IMIifi /2(rp)dr  <  c?l/i"1/2||u||i2(rp)dr.  for  all  v  €  A^p,  i  =  1 j  =  1,...,P7  we  obtain 

A£(/x)  <  ^(/r)cft^1/2Gp(Ai)||u(Ai)  -  (/u)lki©xi,- 

Summing  up  yields 

pr  7 

-  (I/A/’(M)C7t/l"1/2Gp(/7))2^P7||u(/7)  {n)fXi, 

p=  1  i= 1 

and  thus  the  claim.  □ 

5.  Offline-Online  procedure  and  analysis  of  the  computational  costs 

5.1.  PR-SCRBE  approximation 

In  this  subsection  we  briefly  summarize  the  computational  procedure  for  the  PR-SCRBE  approximation 
and  refer  to  [14,  15]  for  further  details. 

Offline.  First,  we  perform  the  pairwise  training  algorithm  as  described  in  [15,  18]  and  apply  a  POD  to  con¬ 
struct  the  empirical  port  modes.  Next,  we  compute  on  each  archetype  component  for  the  source  term  and 
each  coupling  mode  of  the  full  port  space  the  RB  bubble  approximation  by  a  Greedy  algorithm  [26].  We 
remark  that  this  step  is  potentially  rather  expensive  especially  if  we  have  many  degrees  of  freedom  on  the 
ports  as  it  requires  a  couple  of  truth  solves  for  each  RB  approximation  and  hence  each  port  degree  of  freedom. 
Note  that  if  we  allow  ourselves  to  have  only  a  reduced  port  space  at  our  disposal  in  the  online  stage  we  may 
compute  merely  the  RB  bubble  approximations  corresponding  to  the  coupling  modes  of  the  reduced  port 
space  at  possibly  much  smaller  offline  computational  costs.  Independently,  we  have  often  "free"  parameters 
as  Young’s  modulus  or  the  heat  conduction  coefficient  in  which  the  solution  of  (18)  scales  linearly  as  the 
parameters  appear  outside  the  whole  integral .  Furthermore,  we  do  not  have  to  consider  parameters  reflecting 
the  spatial  orientation  as  discussed  in  Section  2.1.  Besides,  the  RB  construction  process  is  embarrassingly 
parallel  and  thus  can  be  performed  in  general  also  for  the  full  port  space  in  a  couple  of  CPU  hours.  Finally, 
we  prepare  the  online  data  sets  and  load  them  to  computer  memory. 


Online.  We  employ  a  graphical  user  interface  to  instantiate  once  I  components  from  the  library,  assign 
the  parameters,  and  connect  the  components  at  ports  to  form  the  global  system.  We  emphasize  that  we 
often  encounter  only  Iejf  <C  7  different  components  in  a  system,  i.e.  components  instantiated  from  a 
different  archetype  or  with  a  different  parameter  configuration,  for  which  we  have  to  perform  the  component- 
local  online  computations.  For  the  remaining  components  we  may  just  copy  the  required  data.  Hence  the 
component-local  RB  approximations  can  be  computed  in  0{Ief  fP1n{QN2  +  N3))  operations,  where  we  omit 
the  superindices  a  and  /  for  Q.  Subsequently  we  assemble  the  component-local  Schur  complement  matrices 
in  G(/e/y(P7)2n2(5iV2)  operations  for  the  Galerkin  formulation  and  0(Ief  f(P1)2n2QN)  operations  for  the 
Petrov-Galerkin  formulation.  Note  that  the  Galerkin  formulation  requires  the  storage  of  the  bilinear  form 
evaluated  in  all  possible  combinations  of  RB  bubble  approximations  for  each  archetype  reference  component, 
which  can  be  critical.  Finally,  we  assemble  the  Schur  complement  system  with  a  direct  stiffness  procedure 
[14,  15]  and  solve  it. 
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5.2.  A  posteriori  error  estimator 

To  realize  an  offline-online  decomposition  of  the  a  posteriori  error  estimator  we  adapt  the  strategy  for 
the  standard  RB  estimator  [7]  and  compute  offline  on  each  archetype  component  the  inner  products  of 
Riesz  representations  (AN (p)  and  AN,1(p))  or  fluxes  (A„(/x))  for  all  (AN (p))  or  some  (AN,1(p)  and  A n(p)) 
combinations  of  coupling  modes,  associated  RB  basis  functions,  and  affine  terms  in  the  bilinear  form.  Hence 
for  A n(p)  we  would  for  instance  compute  the  flux  C\'^ k  by  solving 


L2(t)  ~ 


(78) 


where  •  k  denotes  the  1th  RB  basis  function  and  j  =  1, ...,  P7,  k  =  1, ...,  A fp,  l  <  N,q  =  1, ...,  Qa.  Online 
we  may  then  compute  the  local  indicators  A^ (p) ,  A^ (p) ,  A^p)  by  forming  a  linear  combination  of  the 
precomputed  inner  products.  Alternatively  the  PR  error  estimator  can  be  computed  online  without  offline 
preparations  by  assembling  and  solving  (48).  Next  we  discuss  the  offline  and  online  stages  for  the  a  posteriori 
error  estimator  in  detail. 


Offline.  As  indicated  at  the  beginning  of  this  subsection  we  have  to  precompute  and  store  the  inner  product 
for  the  Riesz  representations  for  all  combinations  of  coupling  modes,  corresponding  RB  approximations,  and 
affine  terms  in  the  bilinear  form  to  accomplish  an  offline-online  decomposition  of  AN(p).  Note  that  the 
necessary  Riesz  representations  have  already  been  computed  (and  stored)  during  the  RB  space  construction 
as  they  are  required  within  the  greedy  algorithm.  Hence  the  additional  offline  costs  are  in  general  smaller 
than  the  costs  for  the  construction  of  the  RB  bubble  approximations.  However,  as  the  required  storage  is 
significantly  larger  than  the  one  required  for  the  data  sets  for  the  Schur  complement  system  the  computation 
of  AN  (p)  is  in  general  especially  feasible  for  systems  of  rather  moderate  RB  space  dimensions  and  moderate 
values  of  Q  and  Mp . 

In  contrast  for  the  computation  of  AN,1(p)  the  data  sets  for  each  Riesz  representative  in  (44)  can  be  stored 
separately.  Moreover,  no  additional  offline  costs  arise  as  the  inner  products  of  the  required  Riesz  representa¬ 
tions  for  the  affine  terms  in  the  bilinear  form  and  the  RB  basis  functions  have  already  been  computed  during 
the  greedy  algorithm. 

To  realize  an  offline-online  decomposition  of  An(p)  we  first  compute  the  fluxes  Clrq^  k  (78)  and  also  the  respec¬ 
tive  fluxes  for  the  source  term.  As  the  right  hand  side  in  equation  (78)  has  already  been  computed  during 
the  data  set  preparation  of  the  Schur  complement  matrix  and  the  left  hand  side  of  (78)  is  restricted  to  the 
ports,  those  additional  offline  computational  costs  are  small  in  comparison  to  the  ones  for  the  PR-SCBRE 
approximation.  To  prepare  the  inner  products  of  the  jumps  of  the  conservative  fluxes  (49)  we  have  to  identify 
groups  of  two  components  which  may  connect  within  a  global  system,  loop  over  each  pair  of  local  ports  within 
each  of  these  groups,  and  form  the  inner  products  of  the  various  fluxes.  We  remark  that  thanks  to  the  port 
compatibility  requirement  and  the  mutual  disjoint  port  assumption  only  a  reduced  number  of  combinations 
of  coupling  modes  and  associated  RB  basis  functions  have  to  be  taken  into  account  and  that  we  hence  expect 
the  required  storage  to  be  comparable  with  the  one  required  for  the  data  sets  for  the  Schur  complement 
system  for  the  Galerkin  formulation. 

Finally,  we  compute  the  constants  c4  and  solve  the  eigenvalue  problem  associated  with  Friedrichs  inequality 
with  negligible  additional  offline  cost. 

Online.  Also  for  the  computation  of  the  a  posteriori  error  estimators  A{p)  and  A1(p)  we  exploit  that 
typical  systems  often  feature  only  Ieff  <C  I  instantiated  components  for  which  we  have  to  perform  the 
component-local  computations.  The  operation  count  for  the  computation  of  AN  (p)  and  A n(p)  is  thus 
0(Ief f{P1)2n2Q2N2).  Although  the  operation  counts  for  A{p)  are  therefore  higher  than  the  ones  for  the 
PR-SCRBE  approximation,  we  expect  the  actual  computational  time  for  the  error  estimator  to  be  within  a 
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reasonable  range  as  we  only  need  to  evaluate  the  sum  over  the  precomputed  inner  products.  To  compute 
AN’1(fi)  we  form  the  inner  products  of  the  Riesz  representations  (44)  in  0(/e//P7?rQ2fV2)  operations  and 
compute  the  sum  in  the  local  indicators  A^v,1(/x)  in  0(Ief f(P’y)2n2)  operations.  For  moderate  values  of 
Q  the  operation  count  for  AN,1(/i)  is  hence  comparable  to  the  one  for  the  computation  of  the  PR-SCRBE 
approximation . 

As  mentioned  at  the  beginning  of  this  subsection  we  may  alternatively  compute  A„(/z)  only  online.  We 
assemble  the  inactive  stiffness  matrices  and  right  hand  side  vectors  in  (48)  associated  with  the  coupling  modes 
ipi,j,k ,  i  =  1,  T,  j  =  1,...,P7,  k  =  n  +  fp  in  0(IeffQNP'rn(Af-p  —  n))  operations.  Thanks  to  the 

^2(7bi)-orth°normality  the  port  modes  the  norm  of  the  jump  can  then  be  computed  in  0(Ieff  P7 n (Afp — n ) ) 
operations.  If  we  consider  port  modes  that  are  not  L2  (ji.j )-ortlionormal,  e.g.  due  to  a  geometric  mapping 
that  transforms  the  port,  we  have  to  additionally  assemble  the  matrix  of  the  left  hand  side  in  (48)  which 
can  be  realized  in  0(P7^j Afp)  operations.  Here,  P7^  denotes  the  number  of  different  port  shapes  within  the 
global  system  either  due  to  a  different  port  type  or  a  different  geometric  mapping.  Thus,  in  this  case  the 
operation  count  for  computing  the  norm  of  the  jump  adds  up  to  0{P2ff Afp  +  /e//P7A/p),  where  cr  depends 
on  the  employed  solver  to  solve  the  linear  system  in  (48).  Unless  we  consider  a  system  with  a  full  port  space 
of  a  great  number  of  degrees  of  freedom,  i.e.  A fp  ~  QP^nN,  the  online-only  strategy  for  the  computation  of 
An(/i)  seems  to  be  preferable  in  both  cases.  Moreover  for  Q{Afp  —  n)  ~  N2,  which  is  reasonable  for  many 
systems,  we  expect  the  operation  count  for  the  online-only  computation  of  An(/z)  for  L2  (7^ )-orthonormal 
port  modes  to  be  comparable  to  the  one  for  the  computation  of  the  PR-SCRBE  approximation. 

Finally,  we  compute  the  minimal  eigenvalue  of  the  Schur  complement  matrix  to  determine  aapp{n)  and  the 
global  constants  c2pp  and  c£app,  if  necessary.  The  computational  costs  are  small  for  instance  compared  to  the 
assembling  of  the  inactive  stiffness  matrix  as  we  may  employ  a  high  tolerance  for  the  solution  of  the  eigenvalue 
problem  and  only  need  to  form  a  sum  of  the  precomputed  local  constants  to  determine  c2pp  and  c* app-  We 
summarize  that  the  online  operation  count  for  the  computation  of  A which  we  have  employed  for  our 
numerical  tests,  is  thus  0(Ieff(P'ynQ2N2  +  (P7)2?!2  +  QNP7n(Afp  —  n))). 


6.  Numerical  Results 

In  this  section  we  verify  both  the  effectivity  of  the  a  posteriori  error  estimator  A1(/x)  and  the  local 
indicators  Av',1(p,),  i  =  1, ...,  /  and  A()(/i),  p  =  1, ...,  Pr,  and  demonstrate  the  applicability  of  the  introduced 
certification  framework  by  analyzing  the  computational  costs.  To  factor  out  the  effect  of  the  geometry  of 
the  archetype  reference  component  on  the  theoretical  results,  we  consider  a  library  of  one  archetype  beam 
component  component  1  as  depicted  in  Fig.  2.  First,  we  consider  heat  conduction,  where  we  focus  on  more 
theoretical  related  issues  as  the  analysis  of  the  effectivity  results.  Next,  we  consider  linear  elasticity  and 
demonstrate  also  for  this  more  demanding  vector-valued  setting  that  our  error  estimator  is  effective  and 
computable  at  small  additional  online  costs.  Our  implementation  is  in  C++  and  based  on  the  finite  element 
library  libMesh  [37,  38].  For  the  solution  of  the  generalized  eigenvalue  problems  we  have  employed  the 
Krylov-Schur  solver  of  the  library  SLEPc  [39] .  For  notational  convenience  we  add  a  subscript  or  superscript 
rel  if  we  refer  to  relative  error  quantities,  i.e.  quantities  divided  by  ||it(/r)||x. 

6.1.  Heat  conduction 

We  consider  equilibrium  heat  conduction  for  the  beam  archetype  component  component  1  (see  Fig.  2) 
with  heat  transfer  from  the  surface  of  the  beam  to  the  ambient  for  a  prescribed  uniform  heat  generation  qdvm ; 
a  combination  of  a  uniform  heat  flux  on  the  bottom  of  the  beam  and  a  volumetric  heat  generation  c/*™ . 
Such  a  heat  conduction  component  might  be  part  of  a  pin  heat  sink  used  for  instance  as  a  processor  cooler. 
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Figure  2:  Archetype  reference  beam  component 

component  1  with  two  ports;  one  at  the  top  and  one  at  Figure  3:  Distribution  of  the  temperature  u  for  system  1 

the  bottom  of  the  beam 


First,  we  nondimensionalize  length  with  respect  to  a  reference  length  wdlm’°  and  we  introduce  a  non- 
dimensional  conductivity  ratio  n  =  kdlm /kdlm’°  for  the  thermal  conductivity  and  a  reference  con¬ 

ductivity  kdlm’°.  Furthermore  we  nondimensionalize  the  heat  transfer  coefficient  hdlm  with  respect  to 
k  Bi/(kdlm’°wdm’0)  for  the  dimensionless  Biot  number  Bi  =  hdmwdlm>° /kdim,  ancj  introduce  the  non- 
dimensional  temperature  u  =  kd%rn'° (Tdlm  —  Td^>ient)/(qdlm(wdlm’0)2).  Finally,  we  introduce  the  non- 
dimensional  heat  flux  qfi  =  and  the  non-dimensional  volumetric  heat  generation  qvoi  =  qd™  / qdlm . 

On  a  physical  domain  f 1  =  (0,  5 L)  x  (—0.5,  0.5)  x  (—0.5,  0.5)  we  may  then  consider  the  non-dimensional  PDE 

—kAu  =  qvoi  over  Q, 
kVm  •  n  +  nBiu  =  0  on  dfl  \  (r^r  U  Try), 
reVw  •  n  =  qf  i  on  Ty, 
u  =  0  on  r  yy  , 


where  L  denotes  the  length  scaling  and  Tjv  and  T yy  the  Neumann  and  the  Dirichlet  boundary,  respectively. 
We  define  the  following  bilinear  and  linear  forms  on  the  archetype  reference  domain  fl  =  (0, 5)  x  (—0.5,  0.5)  x 
(—0.5,  0.5)  for  the  parameter  jl  =  ( k,  Bi,  L,qvoi,qfi ): 


dv  dw 


dv  dw  dv  dw 


(v,w;p,):=(k/L)  ~7^^+kL  —  —  +  —  —  +  BikL 


Iq  dx  i  dx 

f(v,  m)  :=  L  qvoiii  +  /  qfiv, 
J  Q  J  72 


sy  dx  2  dx  2  dx  3  dx  3 


/90\(7iU72) 


where  71  =  {0}  x  (—0.5,  0.5)  x  (—0.5,  0.5)  and  72  =  {5}  x  (—0.5, 0.5)  x  (—0.5,  0.5)  are  the  local  ports  of  the 
component.  We  set  the  parameter  space  V  to  V  :=  [0.5,2]  x  [0.001, 1]  x  [0.6, 1.4]  x  [—1, 1]  x  [—10, 10].  Note 
that  thanks  to  the  choice  of  the  scale  of  the  Biot  number  the  parameter  k  is  a  “free”  parameter  with  respect 
to  the  RB  bubble  approximation  as  it  appears  in  front  of  the  integral  in  the  bilinear  form.  Finally,  we  define 
the  inner  products 


,  „  A .  f  dv  dw  dv  dw  dv  dw 

lv,w)i>:=  /  — ^  — — h  -77 — 77 - 1 — 777 — 777 — , 

jq  dx  1  dx\  dx2  dx 2  dx 3  dx 3 


dv  dw 


dv  dw  dv  dw 


(v,W,  •  d/L)  /_  JA.  +  L  f ^  1 


In  dx  1  dx  1 


q  dx  2  dx  2  dx  3  dx  3 


(79) 


and  consider  the  solution  space  X  =  {u  £  H1{ fl)  :  v  =  0  on  T/y}. 

Unless  otherwise  stated  we  consider  for  this  test  case  the  Galerkin  formulation  of  the  PR-SCRBE  method 
(35).  Thus  thanks  to  Lemma  3.1  we  obtain  well-posedness  of  the  approximation  and  a/'*’ (/i)  >  0.5  =:  aapp(fj1). 
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To  determine  an  approximation  c^pp  for  c*  we  compute  the  largest  eigenvalue  of  (71)  for  L  =  0.5  and  ob¬ 
tain  c*  <  c*pp  «  1.10464  for  all  /i  €  V.  Thanks  to  the  definition  of  the  inner  products  in  (79),  we  have 
c£ (n)  <  1  =:  cxapp ■  To  compute  an  approximation  of  cx ,  we  use  the  simplified  version  of  the  strategy  pro¬ 
posed  in  [34]  and  obtain  cx  <  cf{p)  <  y/(diam(fl)/7r)2  +  (57Xj)2/3  =:  cxpp,  where  L,;,  %  =  1 ,...,/  denotes 
the  length  scaling  parameter  of  the  i-th  instantiated  component  and  cp(/z)  the  constant  in  the  Friedrichs 
inequality. 

Initially  we  consider  a  component  library  library  1  where  we  use  an  underlying  FE  discretization  with 
AT  =  4941  linear  hexahedral  elements  and  a  reference  port  dimension  of  A fp  =  81.  To  generate  the  empirical 
port  modes  we  apply  the  pairwise  training  algorithm1  [15,  18],  compress  the  data  to  6  POD  modes  and  add 
the  constant  mode  to  obtain  7  empirical  port  modes.  For  the  RB  construction  in  the  offline  phase  we  have 
set  relative  tolerance  for  the  Greedy  algorithm  to  10-7. 

First,  we  consider  a  system  1  of  *  =  5  components  with  component  parameters  fi\  =  (0.5,0.01,1.0,0.5,—), 
=  M3  =  M4  =  (1.0,  0.5, 1.0, 0.1, —),  and  /j5  =  (1.0, 1.0, 1.0, 1.0,  5.0).  The  PR-SCRBE  approximation  is 
depicted  in  Fig.  3  and  we  observe  a  rapid  temperature  drop  in  component  5  (in  the  foreground  of  the  picture) 
because  of  the  relatively  high  Biot  number  Bi  =  1.0. 

Analyzing  the  convergence  behavior  of  the  relative  error  ||it(/r)  —  ux {n)\\xl  and  the  relative  error  estimator 
A }.ei{n)  (58)  for  N  =  Nmax  and  an  increasing  number  of  port  modes  n  in  Fig.  4a  we  observe  that  A)eJ(p) 
is  both  a  rather  sharp  and  an  effective  error  bound.  This  also  applies  if  we  consider  n  =  81  port  modes  and 
hence  no  port  reduction  and  increase  the  number  of  RB  basis  functions  N  (cf.  Fig.  4d).  Next  we  consider 
both  a  reduced  number  of  RB  basis  functions  N  and  empirical  port  modes  n.  We  use  either  N  =  8  RB  basis 
functions  and  increase  the  number  of  port  modes  n  (cf.  Fig.  4b)  or  choose  n  =  6  and  increase  N  (cf.  Fig.  4e). 
In  both  cases  A1(/r)  captures  the  behavior  of  the  respective  dominant  model  error  and  is  able  to  detect  the 
error  plateau  due  to  a  too  small  number  of  RB  basis  functions  in  Fig.  4b  or  port  modes  in  Fig.  4e.  The 
effectivity  changes  only  very  slightly  if  the  type  of  the  dominating  model  error  changes.  We  hence  conclude 
that  for  the  considered  test  case  A1(/r)  captures  the  interaction  between  the  RB  and  PR  error  very  well  and 
thus  may  be  employed  to  balance  both  error  contributions  during  the  online  phase,  which  allows  us  to  realize 
significant  computational  savings  as  will  be  demonstrated  below. 

If  we  compare  the  convergence  behavior  of  the  local  RB  indicators  Afr’1(/x,)  (46)  and  the  local  errors 
\W(p)~ux  (m)IU,  for  i  =  1,4,  5  in  Fig.  4f,  we  observe  that  the  local  effectivities  A,fr,1(/q)/|| u(p)  —  ux (/z)||x; 
change  only  very  slightly  for  increasing  N  and  that  for  all  considered  N  the  local  indicators  are  able  to  detect 
the  error  distribution  on  the  components.  Regarding  the  local  PR  indicators  Ap(/z)  (50),  p  =  1,4,5,  we 
observe  that  the  indicators  for  p  =  1  and  p  =  5  capture  the  error  behavior  of  ||zz(/z)  —  ux (/z)||ai©x2  and 
|| u(p)  —  ux (p)\\x5  perfectly  for  increasing  n  (cf.  Fig.  4c).  Due  to  Proposition  4.5  one  might  expect  that 
the  local  indicator  A^/z)  behaves  as  || u(p)  —  ux(p)\\xt !©a'5  for  growing  n,  but  instead  it  shows  the  same 
convergence  behavior  as  \\u(p)  —  ux (/z)||jy4  (cf.  Fig.  4c).  However,  we  emphasize  that  the  observed  behavior 
in  Fig.  4c  is  indeed  preferable  as  the  error  in  component  5  is  mainly  caused  by  the  strong  temperature  drop 
near  the  global  (Neumann)  port  Ts  (cf.  Fig.  2).  To  reduce  the  PR  error  we  should  therefore  increase  the 
number  of  port  modes  at  port  T5,  as  suggested  by  the  local  indicators  A^(/u)  and  A and  not  T4.  We 
remark  that  similarly  also  A2  (/z)  tends  to  behave  like  || u(p)  —  ux {p)\\x3  and  not  like  || u(p)  —  ux (m)||x2©av 
Finally,  we  note  that  the  local  indicators  A([(/z)  are  able  to  predict  the  error  distribution  on  the  ports  well 
(cf.  Fig.  4c). 

Next  we  consider  a  system  2  of  2  components  with  component  parameters  /j-|  =  (1.0,  Bi,  1.0, 0.5,  — )  and 
P2  =  (1-0,  Bi,  1.0,  0.5, 5.0)  and  varying  Biot  number  Bi.  For  n  =  81  and  N  =  10  we  see  in  Fig.  5a  that  the 
local  effectivities  A,fr,1(/Zi)/||zz(p)  —  ux {p)\\xi  change  very  little  if  we  increase  the  Biot  number  and  that  for  all 


1  Following  the  notation  in  [18]  we  have  chosen  Alsamp]es  =  500  and  7  =  2  in  the  pairwise  training  algorithm. 
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Figure  4:  Behavior  of  the  error  estimator  A 1  (f_i)  (58)  and  the  local  PR  error  indicators  A  ()(//)  (50)  on  the  global  port  Tp  and 
the  local  RB  error  indicators  A;V' 1  (//)  (46)  on  component  i  for  system  1;  e  :=  u(fi)  —  (//.), 


considered  values  Bi  £  [0.001, 1]  the  local  indicators  detect  the  correct  error  distribution  on  the  components. 
If  we  consider  n  =  4  and  N  =  Nmax  we  observe  that  the  local  indicator  Af  (//)  associated  with  the  Neumann 
port  captures  the  behavior  of  ||u(/x)  —  u ^  (At)||.Y2  very  well  for  varying  Bi  and  that  again  tends  to 

behave  like  ||it(/z)  —  (//)||xi  rather  than  || u{p)  —  (/x)[|xi®x2  (cf.  Fig.  5b).  A  possible  explanation  for 

the  fact  that  for  smaller  Biot  numbers  we  have  A \{p)  >  A^ (m)  might  be  that  in  this  parameter  range  the 
temperature  drop  is  much  less  localized  than  for  say  Bi  =  0.5  and  hence  more  strongly  affects  the  behavior 
at  port  Ti. 

All  in  all  we  thus  conclude  that  for  this  test  case  the  observed  behavior  of  the  local  indicators  A£,  p  =  1, ...,  Pr 
confirms  the  effectivity  result  of  Proposition  4.5  and  that  the  local  indicators  A)V’1(/ii)  capture  the  local  error 
behavior  very  well.  Moreover  the  local  indicators  are  able  to  predict  the  error  distribution  on  the  ports  or 
the  components,  respectively,  relatively  well. 

To  investigate  the  behavior  of  the  error  estimator  subject  to  varying  L  and  hence  a  non-orthogonal  geometric 
transformation,  we  consider  a  system  3  of  2  components  with  component  parameters  pi  =  (1.0,  0.5,  L ,  0.5,  — ) 
and  /i2  =  (1.0,  0.5,  L,  0.5,  5.0).  Note  that  thanks  to  the  high  Biot  number  the  temperature  near  the  Neu¬ 
mann  boundary  drops  very  quickly  (similar  to  Fig.  2)  and  the  effect  of  the  length  scaling  of  the  beam  on 
the  temperature  behavior  is  therefore  rather  small.  Considering  this  system  thus  allows  us  to  some  extent 
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Figure  5:  Behavior  of  the  local  indicators  (a)  and  A n(/-0  (b)  for  system  2  and  varying  Biot  number;  e  :  =  u(^)~un  (M)- 

(/r) / c^pp  for  system  3  and  varying  length  scaling  parameter  L  (c). 


to  analyze  the  dependency  of  the  error  estimator  and  the  effectivities  on  the  constants  accounting  for  the 
geometric  transformation.  The  results  are  depicted  in  Fig.  5c  and  we  observe  that  for  n  =  4,  N  =  Nmax 
the  effectivity  V1  (^)  / c^pp2 *  slightly  grows  if  we  increase  L.  However  if  we  employ  the  exact  trace  constant 
c4  (Fig.  5c,  square  markers)  the  effectivity  barely  changes,  although  Proposition  4.5  suggests  a  dependency 
on  L.  Setting  n  =  81,  N  =  10  we  see  that  the  effectivity  of  ??1(/r)/c^)p  deteriorates  slightly  more  than 
anticipated  when  departing  from  L  =  1,  but  that  the  deviations  are  still  rather  small.  We  therefore  conclude 
that  at  least  for  the  considered  test  case  A1(p.)  seems  to  be  a  sharp  and  effective  bound  also  for  systems  that 
include  components  subject  to  geometrical  non-orthogonal  transformations. 

Next,  we  analyze  the  dependency  of  r]n(ii)  on  the  port  mesh  size  h.  To  this  end  we  consider  three  different 
libraries  of  component  1  with  linear  hexahedral  FE,  a  fixed  mesh  size  H  =  0.1  in  iq-direction,  and  varying 
port  mesh  size  h  =  0.2, 0.1, 0.05  in  x2-  and  aq-direction  in  reference  coordinates.  We  have  used  c4pp  associated 
with  the  finest  mesh.  For  h  =  0.1  and  h  =  0.05  we  obtained  7  empirical  modes  and  for  h  =  0.2  we  obtained 
5  empirical  modes3.  The  relative  tolerance  for  the  Greedy  algorithm  has  been  set  to  10-7.  Moreover,  we 
consider  the  Petrov-Galerkin  formulation  (39)  for  the  PR-SCRBE  method.  We  consider  a  system  4  of  4  com¬ 
ponents  with  parameter  configurations  Hi  =  fi2  =  M3  =  (1-0,  0.25, 1.0, 1.0,  — )  and  ^4  =  (1.0,  0.5, 1.0, 1.0, 1.0) 
for  which  we  obtain  a  non-constant  temperature  distribution  on  all  ports  thanks  to  the  moderate  Biot  number 
and  the  high  volumetric  heat  source.  The  behavior  of  7yn(/it)  for  h/ 2  divided  by  rjn{[i)  for  h  for  increasing  n 
is  depicted  in  Fig.  6a  and  we  observe  that  the  ratios  are  indeed  bounded  by  the  expected  upper  bound  \f7l. 
As  we  have  observed  smaller  ratios  for  other  systems  as  for  instance  for  system  1,  we  conclude  that  it  seems 
like  the  effectivity  r)n(/j,)  scales  in  h  as  anticipated. 

Finally,  we  address  the  computational  costs  and  consider  to  this  end  again  system  1.  We  note  that  thanks 
to  the  parameter  configuration  /r 4  =  (0.5,0.01,1.0,0.5,—),  iu,2  =  ^3  =  AH  =  (1.0, 0.5, 1.0, 0.1, —),  and 
/x 5  =  (1.0, 1.0, 1.0, 1.0, 5.0)  we  obtain  Ieff  =  3.  First,  we  compare  for  n  =  7  and  increasing  N  the  on¬ 
line  CPU  time  for  An(n),  AJV,1(^t),  and  for  the  PR-SCRBE  approximation,  i.e.  the  computation  of  the  RB 
bubble  approximation,  the  assembly  of  the  Schur  complement  system  and  its  solution.  We  observe  that  the 


2Note  that  we  divide  by  c^pp  to  factor  out  the  effect  of  the  change  of  <ApP  f°r  varying  L. 

,>  We  have  set  the  parameters  for  the  pairwise  training  algorithm  to  Nsarnpies  =  500  and  7  =  2. 
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(a)  dependence  on  h 


(b)  CPU  time  for  PR-SCRBE^1  (/r)  (c)  CPU  time  vs.  error,  A 


Figure  6:  rjn (//.)  for  h/2  divided  by  rjn  ( fi)  for  h  for  system  4  (a).  CPU  time  for  the  PR-SCRBE  approximation  and  1  (//)  for 
system  1  and  n  =  7  and  increasing  N  (b)  and  plot  of  p/.f//)  —  (/dll  v  ’  A*ei(/r)  versus  the  required  CPU  time  for  system  1 
for  n  =  7  and  increasing  N  (c). 


PR-SCRBE  approximation  scales  linearly  in  N  and  that  the  costs  for  the  computation  of  An(/j,)  increase 
only  very  slightly  from  around  0.086  seconds  required  for  N  =  2  (cf.  Fig.  6b).  Moreover  we  see  a  quadratic 
scaling  of  AN,1(n)  which  confirms  the  operation  count  in  §5.  As  the  CPU  time  for  the  computation  of  A1(^i) 
accounts  for  about  75%  of  the  online  costs  we  hence  infer  that  the  proposed  certification  framework  is  indeed 
computationally  feasible  and  that  the  error  estimator  can  be  computed  at  a  reasonable  additional  online  cost. 
The  computational  costs  for  the  corresponding  “truth”  approximation  account  on  average  for  2.25  seconds, 
where  we  have  employed  a  conjugate  gradient  method  with  a  relative  tolerance  of  ICC08  both  for  the  solution 
of  the  Schur  complement  system  and  the  linear  system  of  equations  for  the  “truth”  FE  approximation.  Hence 
for  a  relative  error  ||it(/z)  —  u ^  (n)\\rx  of  10~06  and  an  estimated  relative  error  A4e;(/r)  of  10~°4  we  gain  one 
order  of  magnitude  already  for  this  small  considered  system  (cf.  Fig.  6c).  Greater  computational  savings 
can  be  realized  for  larger  systems  as  considered  in  [18],  especially  for  situations  with  Ieff  <C  I.  We  finally 
note  that  if  we  consider  n  =  7  port  modes  and  prescribe  (say)  a  relative  error  tolerance  of  lO-02  the  error 
estimator  indicates  that  it  is  sufficient  to  consider  N  =  11  RB  basis  functions  and  we  may  thus  reduce  the 
online  costs  to  around  60%  of  the  online  costs  that  arise  for  N  =  Nmax  (cf.  Fig.  6c).  This  demonstrates  that 
we  may  realize  large  computational  online  savings  by  employing  A4(/x)  to  balance  the  error  contributions 
due  to  RB  and  PR  model  reduction. 

6.2.  Linear  Elasticity 

In  this  test  case  we  consider  isotropic  linear  elasticity  for  the  beam  archetype  component  component  1 
(see  Fig.  2).  Again,  we  nondimensionalize  length  with  respect  to  a  reference  length  w*m,°  and  we  nondimen- 
sionalize  Young’s  modulus  Edlm  with  respect  to  a  reference  value  Edlrn’°.  We  consider  a  traction  Tdlm  £  R3 
acting  on  the  tip  of  the  beam,  which  we  nondimensionalize  with  respect  to  Edm’°,  and  homogeneous  Dirich- 
let  boundary  conditions  on  the  bottom  of  the  beam.  We  do  not  take  into  account  gravity.  We  may  then 
introduce  the  dimensionless  displacement  u  =  udlm / vudlrn’° .  Finally,  we  define  the  non-dimensional  tensor  C 
as 

Cijki  =  ^  _  2i /)  2(f_-j-_i/)  1  ^  U  j)  k,  l  <  3, 

where  we  choose  Poisson’s  ratio  v  =  0.3  for  steel. 
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Figure  7:  Magnitude  of  the  dis¬ 
placement  field  u  for  system  5 


N  =  1 

N  =  2 

II 

CO 

Galerkin 

2.942035  •  10"03 

4.908761  •  10”04 

2.310078  •  10"°4 

Petrov-Galerkin 

3.484242  •  lO"03 

-1.736349-  10"U3 

2.222161  ■  10"U4 

N  =  5 

O 

T - 1 

II 

Galerkin 

2.307751  •  10"U4 

2.3077506  •  10"U4 

6.905847  •  10 U4 

Petrov-Galerkin 

2.306615  •  10"U4 

2.307751  •  10"U4 

6.905847  •  10_U4 

Table  1:  system  5:  aapp( Ji)  for  the  Galerkin  (35)  and  the  Petrov-Galerkin  formulation 
(39)  for  increasing  N. 


We  consider  the  archetype  parameter  vector  fi  =  (E,  L,  f2,  r3),  where  the  geometric  parameter  L  allows  for 
a  length  scaling  of  the  archetype  reference  component  domain  f l  =  (0,5)  x  (—0.5,  0.5)  x  (—0.5,  0.5)  in  the 
ii-direction.  Moreover  we  choose  the  parameter  space  V  =  [1.0,10]  x  [0.5,  2.0]  x  [—1,1]  x  [—1,1].  We  may 
then  define  the  following  bilinear  and  linear  forms  on  Q: 


a(w,  v\  fi) 


f(v;  A) 


-([  dw W  dvk  f  dwi  ~  dvk\ 

^  l  /  vp-  C ilkl  /  vp;  x ij kA  0 ..  ] 

\  Jn  dxi  ox i  Jn  oxj  ox i ) 

E  f  dw1  A  dvk  f  dw1  '  dvk 

+  7"  /  W^Cilkll^ - h  EL  /  Ciskt  "XV- 

L  Jn  oxi  Ox  I  Jn  0xs  dxt 


S  7^  M  7^  1, 


where  the  superscript  i  denotes  the  *-th  component  of  the  respective  vector  and  we  assume  summation  on 
repeated  indices.  Moreover,  71  =  {0}  x  (—0.5, 0.5)  x  (—0.5,  0.5)  and  72  =  {5}  x  (—0.5, 0.5)  x  (—0.5,  0.5)  are 
the  local  ports  of  the  component.  Note  that  Young’s  modulus  parameter  is  “free”  as  it  appears  outside  the 
whole  integral  in  the  bilinear  form.  We  also  define  the  inner  products 


f  dv%  dw 1  +  dvl  dw1  +  dvl  dw1 
x  Jn  dx  1  dx  1  dx  2  dx  2  dx  3  dx  3  ’ 


(80) 


(W>«  -E(V£)/a|]K  +  i 


I"  dvl  dw 1  dvl  dut 1 
'n  dx 2  dx 2  da;3  ’ 


and  denote  with  Td  the  Dirichlet  boundary.  Finally,  we  consider  the  solution  space  X  =  {v  £  H1( fl)3  :  v  = 
0  on  TD}  and  obtain  well-posedness  of  the  “truth”  approximation  and  hence  the  Galerkin  formulation  (35) 
of  the  PR-SCRBE  method  thanks  to  Korn’s  inequality. 

In  order  to  determine  approximations  of  the  constants  c4,  cx ,  and  cx  we  proceed  as  in  Section  6.1.  We  choose 
ctapp(^)  as  an  approximation  of  the  smallest  eigenvalue  of  the  Schur  complement  matrix,  where  we  have  set 
the  relative  tolerance  of  the  employed  Krylov-Schur  algorithm  to  0.1.  For  the  discretization  of  component  1 
we  employ  linear  hexahedral  finite  elements  with  J\f  =  3348  degrees  of  freedom  (1116  mesh  nodes)  and  a  ref¬ 
erence  port  space  dimension  of  M-p  =  108.  We  perform  a  pairwise  training  algorithm4  [18]  and  subsequently 


4Following  the  notation  in  [18]  we  have  chosen  Nsarnpies  =  500  and  7  =  1  in  the  pairwise  training  algorithm. 
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apply  a  POD  to  compress  the  so  obtained  functions  on  the  ports  to  24  port  modes.  By  adding  the  six  modes 
that  span  the  space  of  the  rigid  body  modes  we  obtain  30  empirical  port  modes  both  for  the  Galerkin  and 
the  Petrov-Galerkin  formulation.  We  have  set  the  relative  tolerance  for  the  Greedy  algorithm  to  10~7. 
Initially,  we  consider  a  system  5  of  two  components  with  system  parameters  fii  =  (1.0, 1.0, 0.4, 0.6)  and 
/i 2  =  (1.0, 1.0,  — ,  — ),  where  we  remark  that  the  traction  is  prescribed  with  respect  to  reference  coordinates. 
The  Petrov-Galerkin  PR-SCRBE  approximation  for  N  =  Nmax  and  n  =  108  is  depicted  in  Fig.  7  and  the 
bending  in  Xi  and  X3  directions  is  clearly  observable.  In  Tab.  1  we  compare  aapp(/x)  for  the  Galerkin  and 
the  Petrov-Galerkin  formulation  of  the  PR-SCRBE  method  for  n  =  108  and  an  increasing  number  of  RB 
basis  functions  N  with  aA7(/x).  We  observe  that  for  N  =  2  the  constant  aapp(n)  for  the  Petrov-Galerkin 
formulation  is  negative  which  suggests  to  increase  the  number  of  RB  basis  functions.  In  spite  of  that,  we 
see  that  for  N  >  3  both  for  the  Galerkin  and  the  Petrov-Galerkin  formulation  the  approximations  aapp(/j,) 
converge  rapidly  to  a  limit  of  around  2.307751  •  10~04,  which  is  a  very  good  approximation  of  the  coercivity 
constant  aA7(/x).  For  N  =  Nmax  and  n  =  1  we  obtain  aapp(n)  «  2.316429  •  10~04  both  for  the  Galerkin  and 
the  Petrov-Galerkin  formulation  and  the  approximations  merely  change  for  increasing  n.  We  hence  conclude 
that  the  smallest  eigenvalue  of  the  Schur  complement  matrix  may  serve  in  our  certification  framework  both 
as  a  reasonable  approximation  of  (n)  also  for  small  N  and  n  and  an  additional  indicator  for  the  quality 
of  the  RB  approximation. 

Next  we  compare  the  convergence  behavior  of  || u(/i)  —  (/x )\\xl  and  A 4eZ(/x)  for  N  =  Nmax  and  an  increas¬ 

ing  number  of  port  modes  n  for  the  Galerkin  and  the  Petrov-Galerkin  formulation  (cf.  Fig.  8a).  For  both 
cases  we  observe  a  rather  rapid  convergence  of  the  empirical  port  modes,  for  instance  for  n  =  20  we  already 
obtain  a  relative  error  of  less  than  10-4.  Due  to  the  very  small  coercivity  constant  we  obtain  a  rather  large 
effectivity  constant  ??1(/u),  which  has  of  course  to  be  taken  into  account  when  employing  A4(^x)  as  a  stopping 
criterion.  However,  we  see  that  for  both  formulations  the  error  estimator  A  1(n)  provides  an  effective  bound 
as  it  features  the  same  convergence  rate  as  the  relative  error.  If  we  consider  n  =  108  port  modes  and  increase 
the  number  of  RB  basis  functions,  we  see  that  for  the  Galerkin  formulation  the  estimator  increases  if  we 
enhance  N  =  1  to  N  =  2  and  that  the  relative  error  only  drops  slightly  (cf.  Fig.  8b).  Apart  from  that  A4(/x) 
describes  the  behavior  of  the  error  very  well.  For  the  Petrov-Galerkin  formulation  we  detect  a  slight  increase 
of  the  effectivity  ^(/x)  for  increasing  N  (cf.  Fig.  8b).  Considering  n  =  20  port  modes  and  an  increasing  N,  we 
observe  in  Fig.  8c  that  for  the  Galerkin  formulation  and  N  >  3  the  effectivity  is  nearly  constant.  Considering 
the  same  choices  of  n  and  N  for  the  Petrov-Galerkin  formulation  in  Fig.  8d,  we  see  that  ||xx(^x)  —  u %  (/x)llx* 
stagnates  for  N  =  7  and  A 4ei(/it)  for  N  =  6  which  causes  a  slight  increase  in  the  effectivity  rjl{ii).  Apart  from 
that  we  see  that  A4(/i)  captures  the  behavior  of  the  dominant  model  error  well.  Note  that  the  stagnation  of 
both  the  error  and  A4(/x)  in  Fig.  8c  and  Fig.  8d  is  due  to  the  usage  of  only  n  =  20  port  modes  resulting  in  a 
domination  of  the  port  reduction  error.  We  therefore  conclude  that  also  for  this  test  case  the  error  estimator 
A4(/x)  provides  an  effective  bound  which  may  be  employed  to  balance  the  error  contributions  between  the 
PR  and  RB  errors.  The  inferior  sharpness  compared  to  the  previous  test  case  is  due  to  the  small  coercivity 
constant.  Note  that  this  matches  the  theoretical  results  from  the  previous  section  and  is  thus  in  some  sense 
an  anticipated  behavior. 

Comparing  the  approximation  properties  of  the  Galerkin  and  the  Petrov-Galerkin  formulation  for  n  =  108 
(cf.  Fig.  8b)  and  n  =  20  (cf.  Fig.  8d),  we  observe  that  as  expected  the  Galerkin  formulation  yields  a  faster 
convergence  for  increasing  N.  However  by  using  at  most  3  additional  RB  basis  functions  we  obtain  for  the 
Petrov-Galerkin,  which  is  more  favorable  from  a  computational  viewpoint  (see  §5),  the  same  accuracy  as  for 
the  Galerkin  formulation.  We  hence  infer  that  both  formulations  provide  at  least  for  the  considered  test  case 
a  rapid  convergent  approximation. 

To  analyze  the  constants  in  A1(/x)  and  771  (/l*)  due  to  a  geometric  transformation  we  consider  a  system  6 
of  two  components  with  parameter  configurations  /j,  1  =  (1.0,  L ,  0.2,  0.3)  and  /x 2  =  (1.0,  L,  — ,  — )  and  varying 
length  scaling  parameter  L.  We  use  the  Petrov-Galerkin  formulation.  Note  that  in  contrast  to  the  previous 
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(d)  Petrov-Galerkin,  n  =  20  (e)  r)1  (p)/(aapp(p)cA )  (f)  CPU  time  for  PR-SCRBE,  A  1(/z) 


Figure  8:  system  5:  Comparison  of  the  behavior  of  A)eJ  (/i)  and  | ?/(//)  —  ?/A  (//.)  |  for  the  Galerkin  (G)  and  the  Petrov-Galerkin 
(PG)  formulation  of  the  PR-SCRBE  method  (a)-(d);  e  :=  tt(p)  —  v,)" (//);  system  6:  r/1  (//)/(appp(//.)eA )  for  a  varying  length 
scaling  parameter  L  (e);  CPU  time  for  the  PR-SCRBE  approximation  and  A 1  (//.)  for  system  7  and  increasing  n  and  ;V  =  10  (f) 


test  case  the  length  scaling  stronger  affects  the  behavior  of  the  solution.  We  see  in  Fig.  8e  that  for  n  —  20 
and  N  =  Nmax  the  effectivity  ry1(/r)/(aapp(/r)cx)  slightly  increases  for  growing  L.  Employing  c*  instead  of 
c*pp  mitigates  this  effect  but  we  still  observe  an  slight  increase.  Considering  n  =  108  and  N  =  7  we  de¬ 
tect  a  stronger  variation  of  7y1(/r)  for  growing  L.  However  as  the  ratio  between  the  highest  value  of  77*  (/it)  for 
L  =  1.1  and  its  smallest  value  for  L  =  2.0  is  smaller  than  20  the  deviations  are  still  within  a  reasonable  range. 

At  last  we  investigate  the  online  computational  costs  and  consider  for  that  purpose  a  system  7  with  5 
components  and  component  parameters  /ri  =  (2.0, 1.0,  0.3, 0.2)  and  /J.2  =  113  =  ha  =  1^5  =  (2.0, 1.0,  — ,  — ). 
As  in  our  current  implementation  Dirichlet  boundary  conditions  affect  the  component  configuration  we  have 
Ieff  =  3.  Comparing  the  online  CPU  time  of  the  Petrov-Galerkin  PR-SCRBE  approximation  and  A1  (/it)  for 
N  =  10  and  increasing  n  we  see  that  the  PR-SCRBE  approximation  scales  on  average  quadratically  in  n 
(cf.  Fig.8f).  Both  for  A„(/it)  and  AN’1(/i)  we  detect  a  linear  scaling,  which  is  compliant  with  the  operation 
counts  discussed  in  §5.  Thanks  to  the  additional  costs  for  the  eigensolve  the  portion  of  the  online  costs  of 
the  error  estimator  Aw,1(/it)  amounts  to  approximately  82%.  Thus  we  infer  that  the  estimator  is  also  for  this 
test  case  computable  at  reasonable  additional  online  costs. 
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The  computation  of  the  corresponding  “truth”  FE  approximation  required  on  average  14.75  seconds,  where 
we  have  used  a  conjugate  gradient  method  with  relative  tolerance  10~08.  For  the  solution  of  the  Schur 
complement  system  we  have  employed  a  stabilized  bi-conjugate  gradient  method  with  the  same  relative 
tolerance.  The  computation  of  a  certified  PR-SCRBE  approximation  with  a  relative  error  of  8.30531  •  10~06 
took  on  average  0.4  seconds.  Also  for  this  system  we  obtain  a  high  effectivity  ?71(p.)  due  to  aapp(fi)  ss 
1.863954  •  10~05  <  oc^  (n)  ss  1.104764  •  10-04  which  remains  however  constant  for  increasing  n.  Again  the 
tendency  of  the  error  estimator  A4(p.)  to  overestimate  the  error  due  to  the  small  coercivity  constant  affects 
the  applicability  of  A4(/z)  as  a  stopping  criterion  and  hence  limits  to  some  extent  the  computational  savings 
that  may  be  realized  in  the  online  stage. 

7.  Conclusions 

In  this  paper  we  have  introduced  a  new  certification  framework  for  the  port  reduced  static  condensation 
method  (PR-SCRBE)  which  provides  a  rigorous  and  effective  bound.  By  adapting  the  concept  of  conservative 
fluxes  introduced  in  [19]  to  the  PR-SCRBE  setting  we  have  derived  an  a  posteriori  error  estimator  for  port 
reduction  which  is  based  on  local  indicators  associated  with  the  ports.  To  combine  this  estimator  with  the 
estimator  for  the  RB  error,  which  is  based  on  component-local  indicators  that  are  the  norms  of  weighted 
Riesz  representations,  we  have  exploited  the  orthogonality  of  the  component  interior  "bubble"  function  space 
and  the  coupling  modes  on  the  archetype  components.  We  have  derived  upper  bounds  for  the  effectivities 
both  for  the  local  PR  and  RB  indicators  and  the  PR  and  RB  error  estimators. 

The  numerical  experiments  for  heat  conduction  and  isotropic  linear  elasticity  demonstrate  that  the  introduced 
estimator  is  indeed  an  effective  bound  both  for  the  error  caused  by  the  Galerkin  and  the  Petrov-Galerkin 
formulation  of  the  PR-SCRBE  method  and  also  for  systems  containing  components  subject  to  non-orthogonal 
geometric  transformations.  However,  while  the  error  estimator  provides  a  sharp  bound  for  heat  conduction 
problems  the  effectivities  for  the  considered  problem  in  linear  elasticity  are  rather  high  due  to  the  unfavorable 
coercivity  constant.  The  experiments  also  show  that  the  local  indicators  are  effective  and  that  the  effectivities 
change  only  slightly  for  varying  non-geometric  parameters.  Moreover  for  the  considered  test  cases  the  local 
indicators  have  predicted  the  error  distribution  on  the  components  (RB  indicator)  and  ports  (PR  indicator) 
reasonably  well.  This  demonstrates  their  potential  for  being  employed  within  an  adaptive  PR-SCRBE  scheme, 
which  is  subject  of  future  work.  Finally,  the  error  estimator  captures  the  interaction  between  the  PR  and 
RB  errors  very  well  and  thus  allows  us  to  balance  both  error  contributions  in  the  online  stage. 

We  hence  conclude  that  the  introduced  a  posteriori  error  estimator  facilitates  both  an  efficient  error  control 
for  moderate  coercivity  constants  and  an  efficient  and  (locally)  adaptive  choice  of  the  dimension  of  the  RB 
and  PR  spaces.  As  the  numerical  experiments  show  that  the  error  estimator  is  computable  at  reasonable 
additional  online  costs,  the  estimator  thus  enables  us  to  realize  significant  online  computational  savings. 
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Appendix  A.  JT1//2-orthonormal  port  modes 

We  may  alternatively  consider  port  modes  which  are  orthonormal  with  respect  to  a  H1/2- inner  product. 
As  the  computation  of  the  inner  product  corresponding  to  the  Sobolev-Slobodeckij-norm  (cf.  [40])  is  com- 
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putationally  challenging  we  introduce  instead  an  implicitly  defined  H1^2- inner  product,  which  is  defined  via 
an  extension  operator  to  the  archetype  (reference)  component.  In  detail  we  introduce  for  each  local  port  7 j, 
j  =  1,  ...,P7,  an  extension  operator  Vt,  :  A^p  ,  with  Dlyy  =  0  on  7 y  for  j  ^  j' ,  and  then  define  an 

H1/ 2  -inner  product  as 

(C>  '■=  V(,  k  €  .  (A.l) 

The  application  of  the  trace  theorem  yields  that  (•,  -)Hi/2^.yr  is  positive  definite  and  hence  indeed  an  inner 
product  on  the  space  Aj  p .  Note  that  for  the  trace  operator  T,  :  — >  Aj  the  set  {‘Zj&jXj,k}'k=i  is  a 

Riesz  basis  for  AJ^P . 


To  compute  the  respective  port  modes  within  the  pairwise  training  algorithm  we  extract  functions  on 
the  port  for  each  parameter  value  of  the  train  sample5  and  subsequently  compute  the  lifting  to  the  reference 
component.  Note  that  this  step  requires  0(Nsampies[P1  Afp{QN2  +  N3)  +  (P^Af-p^QN  +AfrT})  operations  for 
the  Petrov-Galerkin  and  0(Nsampies[P'yAfv(QN2  +  N3)  +  (P^NvYQN2  +  Afa])  operations  for  the  Galerkin 
formulation  as  we  employ  a  non  port  reduced  SCRBE  approximation  to  compute  the  functions  on  the  ports. 
Here  Nsampies  denotes  the  size  of  the  train  sample  and  a  depends  on  the  employed  solver  to  compute  the 
extension.  We  remark  that  in  general  we  expect  the  computational  costs  for  the  lifting  to  be  rather  small  as 
for  most  snapshots  on  the  port  the  extension  will  have  support  only  within  a  small  part  of  the  component. 
To  guarantee  that  we  obtain  the  full  port  space  we  also  add  extensions  of  Legendre  modes.  The  latter  can  be 
obtained  by  solving  a  singular  eigenvalue  problem  in  the  complement  space  of  the  empirical  modes  [15,  18]. 
Finally,  we  apply  a  POD  to  the  extensions  to  determine  the  principal  components.  Note  that  the  formation 
of  the  Gramian  can  be  performed  in  0(N2amplesAf)  operations  and  that  the  solution  of  the  corresponding 
eigenvalue  problem  requires  0(N3amples )  operations. 

To  define  a  iJ1/2-inner  product  for  a  mapped  port  7 ij  we  recall  our  assumption  that  7”de^  when  applied  to 
a  port  7 j  is  a  pure  rigid  body  motion  which  we  denote  for  notational  convenience  with  7)r5.  We  also  introduce 
the  inner  product  {v,v)xrb  :=  {(J^nav)~1(v  o  Tfj),  (T^nap)~1(v  o  m)x  and  the  mapped  extensions 

ijK  =  (7jnap)(Dljfi  °  (7*yj>)_1)  for  all  k  €  A ^p .  Then  we  obtain 


(iH v  1 


and  may  thus  define 


■Afp 


(C>  K0ff1/2(7M)<ir  (C>  —  i'R/C-  ^jk)x  VC,  n  €  A jj 


(A.2) 


Note  that  thanks  to  this  definition  we  obtain  port  modes  that  are  lV1//2-orthonormal  on  each  global  port  Tp, 
p  =  1, ...,  Pr  within  the  system  in  physical  coordinates.  The  conservative  flux  may  be  defined  as  in  (47),  with 
the  L2-inner  product  replaced  by  the  iV1/2-inner  product  (A.2).  Thanks  to  the  definition  (A.2)  and  the  so 
maintained  orthonormality  of  the  port  modes  the  computation  of  the  conservative  flux  again  reduces  to  the 
calculation  of  the  euclidean  norm  of  the  vector  associated  with  the  residual  in  (48).  If  we  allow  geometric 
mappings  that  deform  the  port  the  computational  costs  increase  as  discussed  in  Section  5.  The  estimate  for 
part  II  in  the  proof  of  Proposition  4.2  can  be  done  in  an  analogous  manner.  Note  that  the  constant  c(  has 
to  be  modified.  The  effectivity  r/n  (//)  can  now  be  bounded  as  follows. 


5Note  that  the  parameter  here  additionally  accounts  for  the  parametrization  of  the  boundary  conditions  on  the  non-shared 
port.  For  further  details  we  refer  to  [15,  18]. 
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Corollary  Appendix  A.l.  For  H1/2 -orthonormal  port  modes  the  local  error  indicators  A/(/i)  may  be 
bounded  as 

A£(/z)  <  l/S ([j,)Cp\\u(fj,)  -  u % (M)l|xi©xi, ,  for  irp  =  {(i,j),  (*',/)}•  (A.3) 

Furthermore  the  effectivity  i)n{p)  satisfies 


t  x 

aapp(ti)  app  app 


Vpf 


max 

p=i,...,Pr 


4- 


(A.4) 


Proof.  The  proof  follows  the  lines  of  the  proof  of  Proposition  4.5  but  does  not  require  the  application  of  the 
inverse  estimate.  □ 
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